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I. INTRODUCTION 

Variational principles are ubiquitous in electromagnetism 0,|2illiSIEIESIlilSEiO|- Thomson's and 
Dirichlet's Principles allow the determination of potentials or capacitances in electrostatics, or the deter- 
mination of steady-state currents and voltages in electrical networks. General reciprocity relations lead to 
variational estimates for eigenmodes of cavities or waveguides, impedances of cavities, apertures, or other 
structures, and even certain scattering coefficients. In magnetohydrodynamics, plasma kinetic theory, 
and certain equilibrium or non-equilibrium thermodynamic situations, stability or dynamics accessibility 
of charged particle distributions has been derived from minimum energy/free ener gy, maximum entropy, 
or minimum entropy production or minimum heat production principles [l^lTsLlMllSllT^lr^lTsLIT^Eoj . 
Density functional theory, an elaboration of the famous Thomas-Fermi and Hartree-Fock methods in 
quantum mechanics, has also been applied to classical and quantum Coulomb gases pH |23 |. Many 
numerical methods and approximation techniques H, |2j 0, |23, mechanical or manual 

computation, associated with the names Raleigh, Ritz, Galerkin, Finite Elements, Minimum Residuals, 
Method of Moments, etc., may all be derived from or interpreted as variational principles. As is well 
known, all of ray optics may be derived from Fermat's Principle of Least Time, and ultimately, all of 
classical electrodynamics may be derived via Hamilton's Principle, a variational formulation demanding 
stationarity of the action functional. 

Variational techniques or approaches enjoy many advantages. They provide unified theoretical treat- 
ments and compact mathematical descriptions of many physical phenomena, which allow for straight- 
forward change of coordinates and incorporation of additional constraints. Conventional equations of 
motion and conservation laws may be derived in an concise manner. Connections between classical and 
quantum limits are often more readily apparent. Variational principles often suggest appealing physical 
interpretations of the behaviors governed by them. They provide a unified, compact framework for 
performing eikonal, ponderomotive, or other averaging techniques. Perhaps most importantly from a 
practical standpoint, they offer starting points for efficient approximation or numerical computation, 
whereby complicated partial differential equations or integro-differential equations may be replaced with 
more tractable mathematical beasts, namely quadratures, ordinary differential equations, systems of 
algebraic or even linear equations, or ordinary function minimization. Although approximate, these sim- 
plified, projected, or parameterized variational solutions may be more easily interpreted or more easily 
computed or offer more insight than the exact forms. 

Here, we derive another principle, which we will call the Maximum-Power Variational Principle (MPVP), 
together with some surrounding mathematical formalism. Although relatively simple in its statement 
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and scope, somewhere between intuitively plausible and obvious, depending upon one's point of view, 
the MPVP may be of some use in the classical theory of radiation, in particular in the analysis of light 
sources relying on spontaneous radiation from relativistic electron beams in undulators or other magnetic 
insertion devices, or more generally, for the approximation of features of other forms of synchrotron 
or "magnetic Bremsstrahlung" radiation. After some suitable generalization, these ideas may also be 
applicable to the cases of Cerenkov, transition, wave-guide/cavity, Smith-PurccU, or other types of 
radiation emitted by classical particle beams traveling through certain media or structures. 

This work arose out of a recent treatment of coherent X-ray generation via a harmonic cascade in an 
electron beam which radiates in the low-gain regime while traveling through a series of undulators. Penn, 
et al. p9ll3Ctl3ll| approximate the modal structure of the radiation in terms of a paraxial mode described 
by certain adjustable parameters. Some of these parameters are subsequently determined directly by 
dynamical considerations, but some remain free, and are determined at the end of the calculation under 
the assumption that the "proper" or "natural" criterion is to choose the values of these parameters which 
maximize the resulting radiated power assuming that mode shape. This procedure seems so physically 
reasonable, appealing, and plausible that intuition suggests that, in some sense, it must be the correct 
approach under the circumstances. However, at least the present authors' intuitions fail to immediately 
suggest a precise derivation or rigorous justification of this power-maximization idea. Then again, it does 
lead us to expect that the governing principles or concepts should be quite fundamental, and applicable 
to more general radiation problems. 

In this paper we pursue these ideas, developing various Hilbert-space approaches to the spontaneous 
radiation from prescribed current sources, culminating in the Maximum-power variational principle 
(MPVP) applicable to such problems under very general assumptions. Although the final results are 
quite simple and intuitive, the mathematics are developed in some detail in a deliberate and pedagogical 
fashion; much can probably be skipped with impunity by all but the most interested readers. 

We first consider the approximate but commonly-encountered case of paraxial fields, where the structure 
and dynamics of the fields and the development of the mathematical results are greatly simplified. 
Motivated by the parallels between the Schrodinger equation in non-relativistic quantum mechanics 
and the paraxial wave equation of classical physical optics, we present an elementary derivation of 
a convenient Hilbert-space formalism for the spontaneous radiation produced by prescribed classical 
sources in the paraxial limit, including a derivation of the MPVP. Guided by these results, we then 
employ Green function and spherical wave treatments of the general three-dimensional radiation fields 
to generalize these results to the case of non-paraxial fields propagating in free space. Specifically, we 
have in mind applications to the spontaneous synchrotron emission from a reasonably well-coUimated 
relativistic electron bunch in a wiggler or similar insertion device, although the results remain valid more 
generally. 

By spontaneous emission, we mean that the trajectories of the charged particles constituting the source 
for the radiation are considered prescribed functions of time, independent of the actual radiation fields 
generated. That is, the particle trajectories are assumed to be determined by initial conditions, externally 
applied wiggler or other guiding fields, and possibly even space-charge effects (quasi-static self-fields, 
either exact including collisions, or in a mean- field/ Vlasov approximation), while any back-action of the 
radiation itself on the particles - either recoil, absorption, or multiple scattering effects - is neglected. 
Thus the MPVP provides an approximate alternative to calculation of the fields via the usual Lienard- 
Weichart potentials or related expressions Q. 

For the case of a wiggler or undulator, this implies that any gain due to ponderomotive feedback and 
dynamic bunching remains small, so the resulting self-consistent stimulated emission component of the 
radiation can be neglected compared to the spontaneous component. The spontaneously-emitted radi- 
ation and its measurable properties (power, angular distribution, coherence, etc.) can then in principle 
be expressed as deterministic functions of the incoming beam phase space profile (including emittance 
effects and any pre-bunching) and the prescribed external fields, without having to solve the equations 
of motion including radiation to obtain completely self-consistent trajectories for the particles. 

Although the radiation is treated classically, this terminology involving spontaneous/stimulated emission 
is standard in the Free Electron Laser (FEL) literature, which, of course, intentionally borrowed it 
from the quantum theory of lasers. The stimulated emission is that which requires (for ponderomotive 
bunching) the presence of both the static wiggler field and additional radiation, either from an external 
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seed (FEL amplifier), from emission earlier in the same wiggler (SASE), or from a previous pass through 
the wiggler in the presenee of a resonant eavity or mirror system (conventional FEL oscillator). The 
spontaneous emission is, by definition, that which occurs in the absence of additional (applied) radiation 
fields, but of course it still requires the static wiggler field. Truly "free" electrons do not radiate because 
they are not accelerated, and in fact energy-momentum considerations would forbid even a solitary 
accelerated electron from radiating if the interaction providing the acceleration did not also allow for a 
momentum exchange with some external matter (ultimately the wiggler magnets or other electrons, in 
the present case). 

Therefore, the description of the wiggler radiation as spontaneous is similar to the atomic case in that 
it occurs in the absence of applied radiation (i.e., real photons), but it differs in that electrons bound 
in atoms, if somehow excited, can then radiate in the absence of any external fields. The static (i.e., 
virtual photon) wiggler field is therefore playing a role analogous to the nuclear Coulomb field. In the 
average rest frame of a siifficiently relativistic beam, the Weizsacker- Williams method of virtual quanta 
may be used, since the Lorentz-transformed static wiggler field begins to resemble an oncoming beam 
of actual photons. Spontaneous radiation occurs when a backwardly-moving, virtual, wiggler photon is 
scattered by an electron into a forward-moving, real photon. Stimulated emission occiirs in a "three- 
body collision" where an electron simultaneously scatters a virtual wiggler photon and real radiation 
photon into two real photons, both in the same mode as the incident photon. 

Throughout this analysis, we also assume that the charge and current densities are not only prescribed, 
but remain localized in space (so that the far-field, or wave zone, is defined) and time (so that certain 
Fourier transforms exist), in a manner more precisely specified in the course of our derivations. We 
could generalize our treatment to additionally include the effects of uniform background conductive or 
dielectric media, but for simplicity we here assume that, apart from its generation by the prescribed mi- 
croscopic sources in a bounded spacetime region, the emitted radiation otherwise propagates in vacuum. 
Generalizations to allow for non-miiform dielectric tensors, representing wave-guides, lenses, windows, 
or other optical devices, might also be possible, but is not pursued here. As we are interested in the 
causal (outgoing) radiation produced by the given sources, we neglect any incident source- free fields or 
fields from other, remote sources. Because we are working in the fully linear spontaneous limit, any such 
fields may be added in superposition at the very end of the calculation. 



II. FUNDAMENTAL EQUATIONS 



Throughout this paper, we denote in boldface type any real or complex vectors, e.g., a or b, and denote 
with carets any real or complex unit vectors such as a which satisfy 

||a||^ (a* -6)^/^ = 1, (1) 
where for any A''-dimensional complex- valued vectors a, b G C^, 

N 

a • 6 = ^ a„6„ (2) 

n=l 

is the usual Euclidean dot product (without any implicit complex conjugation) of the components of a 
and b in some orthogonal basis. Additional notation will be introduced as needed. 

Conforming to widespread, if not universal, convention in beam physics, we work in the (non-covariant) 
Coulomb, transverse, or radiation gauge, where the vector potential A = A± = A{x,t) is chosen to be 
purely solenoidal (i.e., divergenceless), or everywhere transverse in the sense of Helmholtz's Theorem: 

V-A(a;,i)=0. (3) 

Here, a; G M'^ denotes the three-dimensional position and t e M the time, in any one convenient Lorentz 
frame (typically the lab frame, or possibly the average electron beam frame, for wiggler radiation prob- 
lems). In Gaussian units, the scalar potential $ = ^{x,t) then satisfies the instantaneous (i.e., unre- 
tarded) Poisson's equation 



V^$(a;,t) = -4:TTp{x,t), 



(4) 



4 



given the charge density p, while the vector potential the vector potential A satisfies the inhomogeneous 
wave equation 

V2A(a;,i)-^^A(a^,i)--f Ji(a;,t), (5) 

where c is the speed of light in vacuo, and Jj^ is the transverse (solenoidal) component of the full current 
density J, and which, using Q and local charge conservation, can be shown to be given by 

Jx_{x, t) - Jix, t) - J|| {x, t) = J{x, t) - ^ $(a;, t). (6) 

The physical fields are, as usual, related to the potentials by 

E = E^+E^\=~\f^A-V^, (7) 

and 

B = B^ = V X A. (8) 



Note again that we are describing only the fields generated by the prescribed sources p and J . By 
assumption, any externally applied fields (due ultimately to some additional, remote sources) such as 
undulator or bend magnets are not included in the potentials A or 0, but their effects on the trajectories 
of the charge particles constituting the actual radiation sources are assumed to be cither already included 
in J{x,t) and p{x,t), or else are neglected. 

Without some prior expectation for the typical magnitudes of the dominant radiation frequency Wj^d 
or characteristic transverse spot size wq of a radiation beam, it is convenient to re-scale the space-time 
coordinates using the remaining spatio-temporal scales, namely those determined by the speed of light c 
and the plasma frequency ujp associated with the charge density of the known source. That is, choosing 
some reference number density n, ideally characteristic of that of the actual particle bunch, and defining 
the corresponding linear plasma frequency 

were is the electron mass and e the magnitude of the electron charge, we define normalized (di- 
mensionless) variables: a normalized time r = cOpt, normalized longitudinal position ^ = -^z, normal- 
ized transverse coordinates r = x + ry y = ^ {xx + yy) , normalized three-dimensional coordinates 
^ = + r, the normalized vector potential a{C;T) = a{r,£^;T) = ;^^j4(£c,t), normalized scalar po- 
tential (/)(C;t") = = -^<i>{x,t), normalized charge density m(C;t) = l^^ir ■, t) = -^p{x,t), 
and a normalized current density j(C;''") = ji'i'i^i''') = with transverse (solenoidal) part 
J_L = -^J± = 3 - + V_l) ^(j), where V_l = ^ = x-£- + y-£- is the scaled, (geometrically) 

transverse gradient operator, and we also define d = = z-^ + Vj_ as the scaled gradient operator in 
the full three-dimensional space. 

In these variables, Poisson's equation becomes 

r) = [:g, + Vi) 0(r, r) = ^(r, ^ r), (10) 

where = V_l • Vj_ is the scaled (geometrically) transverse Laplacian, and d"^ = d ■ d is the scaled 
three-dimensional Laplacian. The solutions to this equation just consist of the well-known unretarded 
(i.e., instantaneous) Coulomb fields, 

'^(C;r) = -i,/rf3c'i^ (11) 

which fall off like the inverse-square distance from the instantaneous positions of the charges, so they 
do not contribute to the radiation fields, and can be neglected for our purposes. (These unretarded 
Coulomb fields are "attached" to the charges and move with them even under acceleration, so do not 
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include or in any way constitute radiation.) The scaled wave equation for the vector potential may be 
written as 

- a(C;r) = (^^ + Vi - a(r,^;r) = -j^{r,^;r), (12) 

where a also satisfies the gauge condition 

d ■ a(C; r) = (i^ + ■ a{r, r) = 0. (13) 

All of what wc conventionally consider radiation is associated with the fields derived from this vector 
potential, although these fields also include (in the near-field region) non-radiative components needed 
to cancel the non-retarded features of the Coulomb fields. 

To separately analyze each spectral component of the radiation, we perform Fourier transforms in (scaled) 
time, obtaining, for each (scaled) frequency w, the frequency-domain wave equation 

{d^ + a;2) a(r, ^; w) = + Vi + a;^) a(r, ^; uj) = -j^{r, ^; w), (14) 

and frequency-domain gauge condition 

d ■ a(C; w) = (z^ + Vi) • a(r, uj) = 0, (15) 

where 

aiC,uj)^^JdTe'^^a{{C,T) (16) 

and 

MC,UJ) = ^ I dre'^^j^iCr) (17) 

denote the usual Fourier transforms in scaled time, and are in general complex-valued. Because the 
physical, time-domain fields are real- valued, we necessarily have a(C; —to) = a(C; w)*, so we can restrict 
attention to the analytic signal, or positive frequency {u> > 0), components of the radiation fields. We 
assume that the source charge density p and current density j are at least weakly localized in time, to 
the extent that all needed Fourier transforms of the sources and the resulting fields actually exist (at 
least as generalized functions). 

Corresponding to the scaled, frequency-domain vector potential a{<^;ui), we define the scaled solenoidal 
electric field: 

ea ^ SaiC^j) =ioja{C;uj); (18) 

and the scaled frequency-domain magnetic field 

ba = ba{C,c^) ^dx aiCiv). (19) 

The subscript may be dropped when it is clear from which vector potential these fields are derived. For 
later convenience, we also define the scaled, frequency-domain Foynting vector 

Sa(C;'^) -£a(C;^) X 5„(C;cc;)*. (20) 

For any orientable surface S with unit outward normal h, and positive frequency interval < lvq < iv < 
u)i, the quantity 



■p[a](S;a;o, wi) = J duj Ke J dP'an-Sa 



(21) 



may be taken as the scaled, time-averaged (over an optical period) net electromagnetic power, in the 
bandwidth [wqjI^i], which passes through the surface S. 
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In order that the far-field, or radiation zone, even be well-defined, we assume that the physical source 
j(^;r) remains localized in space throughout the relevant emission process, at least in the direction(s) 
of observed radiation propagation. In the general, three-dimensional case, this means that is 
non- negligible only in some neighborhood of some fixed point ^o- However, from 10 and it then 

follows that the solenoidal component J_l(C;''") will not be strictly localized, but must possess decaying 
tails which fall off like 0(1/||C — Cs|P) at sufficiently large distances from ^o- That is, j and j± will not 
both be of compact support simultaneously. But the fields arising from any particular part of j_\_{C', t) 
will consist of the near-zone and intermediate zone fields which will fall off like like 0(1/||C — C'll^) or 
faster, as well as the true radiation fields which fall off more slowly, i.e. as 0(1/||C — C'lD- Because our 
interest resides in the latter, it follows that as long as jXC;"?") is appropriately localized, we can always 
choose some suitably large but finite region beyond which the contributions of the decaying tails of j± 
to the radiation fields may be neglected to any desired level of accuracy, and we can then act as if both 
j{C',T) and jj_(C;''') have compact support. 



By paraxial, we mean that the radiation from the electron bunch in the insertion device can be assumed 
to be highly mono-directional, in the sense that the relevant angular scales - namely any possible overall 
angular offset between the electron beam and optical axis of the insertion device, the angular spread of 
the particle beam, and the characteristic angle of diffraction - are assumed to be small, in a sense made 
precise below. If we restrict attention to radiation observed within a sufficiently small detection solid- 
angle, then radiation from relativistic particles in certain non-straight devices, i.e., synchrotron radiation 
in circular rings, may also be treated within the paraxial approximation. The needed local directionality 
is provided by the assumed finite acceptance entendue of the detector and the characteristic 0(1/7) 
angular spread in the radiation resulting from relativistic 0(1/7^) compression effects between emitter 
and observer times, where 7 is the usual relativistic kinematic factor associated with the mean velocity 
of the source charges, which together limit the contributing portion of the electron trajectory to some 
small arc. 



We have already seen how, in the Coulomb gauge, the scalar potential (p does not contribute to the 
radiation fields. The vector potential a contributes to near-field (quasi-static), intermediate (induction- 
zone), and far-field (radiation-zone) effects (and in fact part of it must cancel, in the resulting fields, the 
non-retarded effects of the scalar potential), but considerable simplification is possible if we effectively 
restrict our attention to those portions of the fields which, if allowed to evolve freely after propagating 
beyond the finite support of the sources, actually constitute radiation, and further assume that the 
radiation of interest propagates nearly along the +z axis, with a small characteristic angle of diffraction 
9o -€1 1 and any overall angular offset SOq < 0{9d)- Then we can make the so-called paraxial approxima- 
tion to the wave-equation for the vector potential a, which in effect consists of an asymptotic expansion 
in S9o- 

First, for any w > of interest, we express the vector potential in enveloped form: 



where the (scaled), axial carrier wavenumber k = k{uj) = +uj satisfies the vacuum dispersion relation for 
a local plane wave traveling along +z. Similarly, we define an enveloped source term such that: 



III. PARAXIAL CASE 



A. Paraxial Approximation to the Wave Equation 



(22) 



j^(r,^;c.) = 2fc/(r,C;^)e'' 



(23) 



Using the chain rule, the wave equation then becomes 




(24) 
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while the radiation envelope ^p and source envelope / both satisfy envelope-transversality conditions: 

(i^+zfci + Vi) •t/>(r,e;t^) = 0, (25a) 
(z^+zfci + Vi) •/(r,^;a;) = 0. (25b) 

We assume that the charge density p and full current density j remain weakly localized transversely, in 
the sense that they are square-integrable in any transverse plane. Actually we will impose the slightly 
stronger constraint that |/(r;^;w)| more rapidly than l/|r| , as |r| — > oo. We also assume that 
the physical currents are strongly localized longitudinally, in the sense that they are negligible outside 
some finite interval ■^q — ^ — ^-long the propagation axis, implying in general that the solenoidal 
part of the current density j± will actually have a component with non-compact support, corresponding 

to the — (^z-^ + ■§^4' term. However, this contribution to the sources falls off no slower than 

the inverse-square of the distance from the source charges, so with arbitrarily small relative error in 
the radiation fields, which are linear in the sources and fall off with the inverse distance from them, 
we can safely truncate the source envelope /(r;^;u;) outside of some sufficiently large but finite range 
Co < 'Co — ? ^ < 'Ci- (Again, the neglected contributions are strictly necessary to cancel the unretarded 
part of the Coulomb fields to ensure that the full fields remain compatible with relativistic causality, but 
we are confining attention only to actual radiation.) Most of the final results will not actually depend on 
the size of this interval, so if necessary we may safely take |Ci — Col ^ oo at the end of the calculation, 
provided the charge and current densities continue to fall off sufficiently rapidly to ensure convergence 
of the relevant integrals, and we can somehow distinguish the far-fields from the local quasi-static fields 
when the latter may then extend to infinity. 

The (scaled) waist of the radiation beam, or characteristic transverse spot size wq, may be defined in 
terms of the minimal, power-weighted, root mean-square radius: 

T^Wn = ml — 5 — . (26) 

The corresponding characteristic diffraction angle 9^ may then be taken as 

^o = ^- = r^, (27) 

27r kwq 

and finally the (scaled) Rayleigh range Crj or longitudinal length-scale for appreciable change in the 
transverse spot size due to diffraction, becomes: 

Ch = 5^«^o- (28) 
We now treat 9^, as a small parameter. Specifically, if ^pj is any (non-zero) component of the radiation 



/m - O(C-i), while |lV^^,|i/|^,| 



envelope ip, then we may consistently assume that: 

0{w^^). If the spot size is sufficiently large compared to the carrier wavelength A = 27r//c, such that 
^ 1; then with a relative error of 0{9^), we may drop the higher-order longitudinal derivative in 
equation (|24|l . obtaining the forced paraxial wave equation 

0^ + Jk^l) ^(^.C;t^) + /(r,C;t^) = 0. (29) 



Note that we have not yet directly imposed the same ordering assumptions on the spatial variations in the 
envelope f{r; C; i^), but as the latter is the actual source for the radiation, it is clear that the transverse 
scale-lengths of "0 must be consistent with those of /(r; C; w) under diffractive propagation (and therefore 
typically comparable in the vicinity of the sources.) However, as long as we consider sufficiently small 
wavelengths, the scale-length for longitudinal variations in f{r;£^;uj) need not be comparable to Cr 
in order for the paraxial approximation to remain accurate for the fields of interest, i.e., for radiation 
observed in the far-field region beyond the sources. Specifically, we may imagine the paraxial fields as the 
superposition of those fields from each infinitesimal transverse slice of current, evolved forward according 
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to the free-space paraxial propagation, so that each component will satisfy the paraxial conditions 
everywhere except for a discontinuous jump at its source slice, which can be replaced with an equivalent 
boundary condition applied to the homogeneous (source-free) equation. The paraxial solution is expected 
to remain accurate for any ^ > if the paraxial conditions both hold for the calculated fields in the 
region ^ > and also hold for ^ < ^i, when the fields are extrapolated backwards via the free-space 
(homogeneous) paraxial propagation. 

This paraxial wave equation is valid when the propagation direction is sufficiently close to the fiducial 
(-1-^) axis, and the scale-length for transverse variations in the radiation beam are large compared to 
the carrier wavelength and small compared to the length-scale for longitudinal envelope variation (i.e., 
residual variation after the faster carrier oscillations have been removed.) Actually, although we have 
formally assumed that ^ 1, numerical solutions to the full wave equation typically reveal that the 
paraxial approximation actually remains surprisingly accurate for focused radiation beams provided only 
9d < 0{1). For still smaller spot sizes (tighter focus), the paraxial equation predicts a coherence volume 
Vcah < 0{X^) for a single oscillatory mode, which would be in conflict with the diffraction limit dictated 
by the "Heisenberg" inequality constraining conjugate Fourier transform pairs. 

From our above ordering assumptions, and the transversality condition H25a|l . it follows that {ipzl /{\'4'x\^ + 
I ^2/1 y^"^ ~ 0{9j^), so the axial component ipz can in certain circumstances be neglected outright, but 
we will retain it here in order to explicitly satisfy the gauge condition. However, whenever convenient, 
the gauge constraint for fully paraxial fields may be simplified, within the same 0{el) relative accuracy 
as the paraxial wave equation itself, by dropping the longitudinal derivative: 

{tkz + Vi_)-ip{r,C,uj)^0. (30) 

So if the (geometrically) transverse components xp± = rp — zz-xp are specified, then in this approximation 
V'z(r, ^; w) may be taken as 

z-xPir,^;uj) = j:V^-'iP^{r,^;Lo), (31) 

or conversely, if if^zir, lu) is given, then a ■4'± can always be chosen to satisfy (pTT|l and ensure paraxial 
Helmholtz-transversality. The solenoidal source envelope / must in general satisfy the complete envelope 
constraint l|25b|l including the longitudinal derivative, because its ^-dependence may be arbitrary. 



B. Hilbert Space Formalism 



Mathematically, the homogeneous (i.e., free-space, or source- free) part of H29(l is exactly analogous to 
the Schrodinger equation, in the position representation, for a quantum mechanical, spin-one particle 
of mass M, moving non-relativistically in a potential-free, two-dimensional, Euclidean space, where 
the longitudinal position ^ plays the role of the "temporal" evolution variable, r serves as the spatial 
coordinates, and polarization is analogous to spin angular momentum, all in units where h — I and 
M = k = uj. 

In order to leverage the analogies to ordinary quantum mechanics, we introduce a Hilbert space, inner 
product, state vectors, operators, etc., using a Dirac-like notation. The full Hilbert space may be 
considered a tensor product of the Hilbert spaces for the (geometrically) transverse-spatial and the spin 
(polarization) degrees of freedom. (Our notation is similar, but not identical, to that in 32, 33] ). We 
associate the ket (really a spinor) with the complex vector field ipif) defined on the (geometrically) 
transverse spatial plane, i.e., as a mapping ■0 : — > (possibly also parameterized by additional 
continuous variables such as ^ and to and by discrete mode indices, all suppressed for the moment). The 
corresponding bra (■01 is associated with the dual, or conjugate transpose, vector field, xlj{r)\ and a 
combined £<2/Euclidean inner product is defined by 

(-01IV2) = j d^r Mr)* -Mr)- (32) 

We define a number of Hermitian operators acting on the spatial degrees of freedom by their action on 
these position-representation wave vector fields xp. Given any fixed unit vectors e, e' g in the complex 
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linear span of x and y, the transverse position operators just act via multiplication by the corresponding 
transverse spatial coordinates 

e-QV(^) = (e-r)i/j(r), (33) 

while the conjugate "momentum" operators act by differentiation: 

e' ■P^{r)=l{e' ■V^)^{r), (34) 

so these operators then satisfy the usual canonical commutation relations 

[e-Q, e' ■ P]=i{e-e')Ir, (35) 

where Ir is the identity operator on transverse position-state space. We stress that these operators act 
only on the transverse spatial coordinates, so that z-Q = z- P = 0. The "kinetic energy" or "free-particle 
Hamiltonian" operator is then defined as: 

H^mP-P = -l-k'^l' (36) 

which is just the infinitesimal generator of free-space, diffractivc propagation along ^. With somewhat 
more mathematical care about operator domains than is justified for this presentation, these Hermitian 
operators may be extended to fully self-adjoint operators without difficulty. 

In order to include the polarization, we also define the complex helicity basis in : 

is = -s^ (x - isy) + (1 - s^)z = -^5^(1 + s)e^ + ^^{1 - s)e+ + (1 - s^)z, (37) 
for s = — 1, 0, +1, satisfying 

e* • is' =Sss', (38) 

where Sg s' is the usual Kronecker delta symbol and e± = -^(i ± iy) = e'^, and then can define various 
Hermitian "spin" operators by their action on the polarization components of the fields ip. In particular, 
we define a vector of mutually orthogonal polarization projection operators 

n= J2 (39) 

s=-l,0,+l 

where 

n,* = e, (e: • *) , (40) 

such that 

n^Hy = n,/n, = 5s«'n„ (4i) 

and 

E n. = /3, (42) 

s=-l,0, + l 

where 

s=-l,0, + l 

is the identity operator on the spin degrees of freedom. We also define the helicity, or circular polarization 
operator by 

5, = ni-n_i, (44) 

so that 

S^ip = e-e*_-ip-e+el-ip. (45) 
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and 

Sl = h-Iio. (46) 

For an enlighte ning discussion of the orbital and spin components of angular momentum in electromag- 
netic fields, see |3^. 

Clearly, all of these polarization operators commute with the spatial operators, so we may define gen- 
eralized (i.e., non-normalizable) simultaneous eigenkets of position and spin |r;s) for any r e and 
se{-l,0,l}: 

Q\r-s)=r\r-s), (47) 



e''--^|r;s) = |r-r';s), (48) 

and 

S,\r;s)^s\r-s). (49) 

We similarly define generalized simultaneous eigenkets of spin and "momentum" (transverse wavevector) 
by taking the Fourier transforms of these position kets: 

\p\s)^^ j d^re^P-^\r-s), (50) 

which satisfy 

P\p-s)^p\p-s), (51) 

and 

S,\p-3)=s\p-s). (52) 
These generalized eigenkets satisfy the orthogonality relations 

{r-s\r'-s')^5,s'5^^\r~r'), (53) 

and 

{p-s\p'-s')^5ss-5(^\p-p'), (54) 

and the completeness relation 

^ j cfr \r-s){r-s\^Y. j d'p \p; s){p; s\ ^ I , (55) 

where (r) is the two-dimensional Dirac delta function, and / = /s (g) is the identity operator on the 
full Hilbert space H consisting of all which are normalizable, i.e., for which < oo, implying 

that the corresponding paraxial radiation fields '4>{r) are of finite power (but not of finite energy). The 
allowed radiation envelopes actually reside in the Hilbert sub-space H± C H consisting of vector fields 
which also satisfy the transverse envelope gauge constraint when evolved in ^ according to free-space 
paraxial evolution. 

Defining the coupled spin-momentum operator 

T{5) =^ (k~i6-^)llo + P-Il (56) 



for a real parameter 6, and the constant reference ket \tpo) G ^ such that 
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the full envelope-transversality constraints (|25|l on the radiation and source may be written in our 
quantum notation as: 

{^Po\TiS^^)\^P{C;Lu))^0, (58a) 

{iPo\T{S^l)\f{(;cj))^0, (58b) 

while the approximate paraxial gauge constraint H3()(l becomes: 

{^Po\T{S^0)\1P{tu;))=0, (59) 

Again, within the paraxial approximation, either of these constraints may be used interchangeably for 
the radiation, but the source must satisfy the full transversality constraint because its ^ dependence is 
arbitrary. 

It is straightforward to see how the full envelope gauge constraint (|58a|l (including the longitudinal 
derivative) will be automatically maintained by the evolution, and how the approximate paraxial con- 
straint H59() (neglecting the longitudinal derivative) will be preserved by free-space propagation or by 
propagation in the region of a fully paraxial source. Assuming that \^p{£^^,Lu)) satisfies the paraxial wave 
equation, with an initial condition satisfying 

hm (^Po\T{S)\^P{(■,iu))^0, (60) 

(our specific choice of \ip{£^;Lu)) — for all ^ < Co clearly works), and that the source \f{^;Lj)) satisfies 

(V'o|T(<5)|/(^;c.)) =0 (61) 

for all ^, and noting that 

[T{5),H{i.)]^\T{6),^]=0, (62) 



^ff (^) = 0, (63) 

and 

HiLo) \xPo) = 0, (64) 

then: 



= -I {rPo \mH{uj)\i,{e, oj)) + I {^0 \T{S)\ /(^; uj)) (65) 
= -l{^Po\Hicu)T{S)\^P{(■,tJ))+0^0. 

Using the initial condition, we then have {xpo \T{6)\ V'(Ci ^)) ~ for all ^. This shows that the full {6 — 1) 
envelope transversality constraint (|58al (including the longitudinal derivative) on the fields is preserved 
everywhere, if it holds initially for the fields and holds everywhere for the sources. The approximate 
(i.e., 6 = 0) paraxial gauge constraint H59|l is preserved by the evolution in regions free of sources and 
wherever the solenoidal source also satisfies this same, stronger condition. 



C. Green Function Solution 

In this quantum-like notation, the paraxial wave equation may be written as 

^ Hm^Lo)) -imu;)) (66) 
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Given some "initial" state , the solution to the homogeneous part of the equation (i.e., for 

|/(^; Lo)) = 0) representing free space propagation of the radiation fields, may be written in terms of the 
unitary evolution operator as 

m;u;))^U{i,^'-cj)m';u;)), (67) 
where U{^,S,'; k) satisfies the operator- valued Schrodinger equation 

ij-^U{^,^'-u)^HU{i,i'-uj) (68) 

with initial condition 

C/(C',r;a;) = /, (69) 
and possesses the group composition properties 

(7(e, e'; u)-^ = [/(C, e'; t^)^ - U{i', u^), (70) 

and 

u{^,e';co)U{e',e;^) = u{^,^';io). (71) 

Here H = H{u!) (with uj = k) is the Hamiltonian, or diffraction operator, defined above. Because H is 
^-independent, we can immediately integrate (|68|l to find 

[/(e,e';c^) =e-*«-«')^("). (72) 

The solution to the full inhomogeneous wave equation (|66(l may be written in terms of the propagator, 
or causal Green function operator G(^,^';w) satisfying 

(z^-i/)G(e,e';c.) = <5(e-n (73) 

and 

G(C,^';a;) = for ^ < (74) 
It can easily be verified that this operator is given by 

Git e-, = - ouit e-, (75) 

Here is the usual Heaviside step function, and S{£,) ~ ^©(■C) is the one-dimensional Dirac delta 
function. Assuming no initial (incoming or outgoing) radiation, the formal solution to the inhomogeneous 
equation (i.e., with the source term) is then 

00 

m;u^)) = -[ dC'G(C,C';^)l/(^';^)) 



min(^,^i) ("76) 

= 1 I <[/(C,C';c.)|/(^';a;)), 

where we have used the assumption that the support of /(^;r;w) is confined to < ^ < Ci- For smy 
i < Co, we then have \if>{t,uj)) = consistent with our assumed initial condition, while for all ^> £,1, this 
just reduces to free-space propagation: \^p{£_■,uJ)) = U{£_,£,i;lj) \ip{£i;uj)) . In between we must actually 
solve the full expression lf7B|) . 
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D. Energetics 



In our dimensionless variables, recall that the normalized, solenoidal electric field is given by 

£a±(r;C;T) = -|ra(r;^;r) = -|:V,(r; ^; r)e*'=« (77) 

in the (scaled) time domain, or equivalently 

Sa±{r;^;ui) = iu a{r;^;ui) = iuiip{r;^;u))e''''^ (78) 

in the (scaled) frequency domain. Within the paraxial approximation, the inner product 
{ip{^;oj)\ip{^;uj)) is then proportional to the (time-averaged) power spectral density for radiation of 
frequency lv passing through the transverse plane at the longitudinal position ^. In fact, in scaled units 
we may define the (normalized) power spectral density as 



^T,M[a](^;a;)=a;2(TA(C;a;)|T/.(^;a;))=a;2y' |V(r;^;a;)f 
= w^J (fr \a{r;^;w)\' = J d'r \e^{r;^;iv)\' 



(79) 



Strictly speaking, this yields the time- averaged Poynting flux only to O(^^d) at any finite distance from 
the sources, but it is exact as ^ ^ oo and we consider only the truly radiative transport of energy to 
spatial infinity. The unitary (inner-product preserving) nature of the paraxial evolution in vacmun then 
corresponds to energy conservation for paraxial fields, where the (time-averaged) power crossing any 
transverse plane in any frequency bandwidth, remains constant (with respect to propagation distance 
^) during free-space evolution. Note, however, that while it is particularly convenient to formulate 
paraxial propagation in terms of the transformation of the fields in successive transverse planes, the 
actual wavefronts (level surfaces of phase) of a general paraxial beam are not planar except at the focus. 

Prom the actual solution \tp{£^2',i^)) in any conveniently-chosen transverse plane ^2 > ^1 beyond (i.e., to 
the right of) the sources, we can construct the extrapolated free-space field 

|x(e; uj)) = U{i, 6; a;) |t/;(6; u^))=ij < U{i, e'; a;) m'; a;)) , (80) 

io 

which represents the post-source paraxial radiation envelope extrapolated throughout all space via free- 
space propagation, including in the backward or "upwind" direction into the region where sources are 
actually present. Thus \x{^\ '^)) and |V'(^; t^)) coincide for all ^ > ^1, but not in general for ^ < ^1, where 
|x(^; uj)) satisfies the homogeneous (source- free) Schrodinger equation with no incoming radiation as ^ — > 
—00, whereas \'ip{S^;uj)) satisfies the inhomogeneous equation everywhere, so that all outgoing radiation 
at ^ ^ 00 is matched by incoming radiation at ^ —00. For any particular ^, these extrapolated fields 
are those fields which would have been present if the fields actuajly observed beyond the sources at some 
^2 > were instead produced by some effective source 9(t";^;^;w) located in some sufficiently remote 
region ^0 — ^ < C < — ^, instead of by the actual sources within ^0 < ^ < : 

€i-f 

\x{^;uj))=i J d^' U{^,^';u;)\g{^';^;u;)), (81) 
£o-f 

where the effective source is determined from the actual source simply by longitudinal translation and 
compensation for diffraction: 

\g{^'; e, u;)) = U{^', ^ + |/(^' + co)) , (82) 
and ^ = ^(^;Ci) > may any positive constant satisfying 

^>ei-e (83) 
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Choosing any ^ > ^1—^07 the expression H81|l is then valid for any ^ > or we may actually take 
^ —^ +00 if convenient, moving the effective sources to an abitrarily remote "upwind" location. 

Using the Green function solution H76(l . and taking advantage of the assumed finite support of /, we 
find that for any ^ > $1, 

min({,{i) 



(0 



W {x{e;^)\f{e;u;)) ^^^^ 



which is just a three-dimensional inner product, or overlap integral, between the source envelope 
/(r;^';w) and the free-space radiation envelope x(r;^';w) which has been extrapolated via free-space 
propagation backwards into the region where sources are present. 

While solenoidal and irrotational vector fields are not in general locally orthogonal in the geometric 
sense, they are Hilbert-space orthogonal when their Euclidean inner product is integrated over all space. 
Because \xp{£_;uj)) and hence \x{d^)) satisfy the transverse gauge constraint, this overlap integral may 
be simplified to 



{^Pi^■,L,)\^P{C;Lo)) / d^J d'r x(r-; C'; a;)* • /(r; ^) 

^ J d^J d'r [^u;xe^''T ■[2kfe^^^] 



J dej d^r e^^ir;^;^)* ■j±{r;e;uj) (85) 
J d^'J d'r £^^(r;e»*-j(r;^';^) 

/ de (£^j.(r;^)U(r;^)>, 



where e-^±{r; ^' ;lu) is the scaled, solenoidal electric field associated with the extrapolated radiation 
envelope X- In order to evaluate the power, we therefore need never explicitly decompose the current 
density into its solenoidal and irrotational components, nor worry about the precise cutoff beyond which 
the solenoidal component may be neglected. 

As one would expect from energy conservation considerations, this overlap integral may be interpreted 
in terms of the rate of energy transfer between the charges constituting the sources and the fields. We 
may define the scaled power spectral density associated with the rate of (time-averaged) mechanical 
work done by scaled electric fields e(r;f;w) on the charges associated with the scaled current density 
j{r;^';uj) as 



j] (c) = Re I < {s{e;uj)\j{e;io)) 

= ReJ di'j d^r e{r;^';u:)* -jir-i'-Lo). 



(86) 



From (|85|l . we see that j d£/ {£~^^{£^' \Ld)\j{£^' must be purely real and non-positive, because {ip\ip) = 
II ■011^ is necessarily real and non- negative. Physically, this reflects the fact that the back-extrapolated 
fields, if actually present, would lead to no reactive energy transfer to local fields associated with re- 
arrangement of the source charges, and mechanical energy would be transferred monodirectionally from 
the charges to the radiation fields. Therefore H85|) is equivalent to: 

^T,M[a](e;c.) = -i^jT^Je^^; j](c.) (87) 
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That is, the average power seen in the far field is exactly one-half of that which would be delivered by the 
sources to the extrapolated homogeneous fields if these fields were those actually present in the region 
of the sources. 

The presence of the extra factor of ^ is intuitively plausible. The actual paraxial fields from each 

transverse current slice consist only of the forward half of the corresponding extrapolated fields from this 
slice, and hence the former fields could interact only with "downstream" current slices while propagating 
in the forward (+z) direction, while the extrapolated fields from any current slice impinge on all other 
current slices, both upstream and downstream. If we adopt the symmetric convention 6(0) = 5 for 
the Heaviside step function (which is anyway the natural choice when Fourier transforms are involved), 
each ciirrent slice also appears to interact only with one-half of its own field. Assuming a reciprocity 
arising from Newton's Third Law (magnetic forces can violate this law, but do no mechanical work, and 
radiation reaction forces can violate the third law as well, but are ignored here), we then anticipate an 
over-counting of the work by exactly a factor of 2 when we use the extrapolated fields and integrate over 
all current slices. This is a dynamical, electromagnetic analog of the well known electrostatic result that 
the potential energy associated with a given charge distribution p in a given external potential $ext is 
given hy U = J cPx /5$ext, while the self-energy (potential energy of formation) of a charge distribution 
resulting in a potential ^ is U = ^ J cfix p<I>. 

Below, we will provide an elementary proof of this result in the general case, but it is informative to verify 
it with an explicit calculation for paraxial fields. By using the Green function solution and inserting a 
resolution of the identity, we have, for any ^ > ^1 , 



_§_ 



€0 io 



.-J: d-r 



(88) 



de {r;s\U{^o,e;^)\f{e;<^)) 



^y^V J d^'g{r;^o;e-^o;uj] 



which explicitly demonstrates its reality and non-negativity. Using the same sort of manipulations, we 
then have 



= Re 



So 

-2iivk J di' (t/'(0)|/(€';a;)) 



«0 

<|< {m";^)\U{(.",(.';u;)\f{^';u;)) 

^0 «0 



(89) 



r «1 « 

= -2a;2 ^ j d^r Re J d^ g{r; ^o; - 6; oj) J d^" g{r; Co; T - ^0; oj)* 
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But using integration by parts, this simplifies to 

_d_ 



So, indeed it is the case that: 



(90) 



(91) 



consistent with expectations of energy conservation. Note that while the extrapolated source-free solution 
X coincides ,by construction, with a in the downstream region £, > ^i, and so exhibits identical out- 
flowing power, in order to actually satisfy the source-free equation, it must have an equal magnitude of 
power flowing into the interaction region through any upstream plane at ^ < which accounts for the 
extra factor of two. 



E. Basis-Set Approximation 



As we have seen, equation H76|l in principle gives the exact solution to the paraxial wave equation with 
prescribed sources, without any further approximations, and is compatible with energy conservation for 
paraxial electromagnetic flelds. For given sources, the field at an arbitrary longitudinal plane ^ can 
be determined via a convolution integral over the longitudinal source position ^' and a two-dimensional 
integral over either the transverse position-dependence (real space) or transverse wavevector-dependence 
(reciprocal space) of the solenoidal component of the current. In position space, we have 

m;u;))^Y. J d'"- J d'r' J (r'; .|/($^ u;)) \r;s) , (92) 

where we have used the well-known result for the two-dimensional, position-space, free-particle propa- 
gator in quantum mechanics. In momentum space representation, the propagator is diagonal, and we 
have 

min({,{i) 

m;u;))=zJ2jd^p J dC'e-'^(«-«')^(p;s|/(e';^)) (93) 

where (p; s\f{^'; uj)) is just the Fourier transform (with respect to the scaled transverse wavevector) of the 
solenoidal current density slice at longitudinal position Now, using either either of these approaches, 
we require something like 10 integrations to determine a given frequency component of the field at a given 
spatial location, starting from the full current density as a function of space and time coordinates. In 
practice, we often only know the source probabilistically, so we will typically have to perform additional 
averages over the particle distribution function, involving up to 6 more integrations over the full particle 
phase space. 

A complementary approach is to decompose the paraxial radiation fields into a complete, orthogonal 
set of modes. Because we are actually only interested in the radiation observed beyond the region of 
interaction with the localized sources, these modes are naturally chosen to satisfy the homogeneous 
(source-free) paraxial wave equation for all longitudinal positions ^ > ^i, but we may then also choose 
these modes to be extrapolated free-space solutions everywhere, including within the support of the 
sources. That is, we consider a countable set of explicitly ^-dependent, solenoidal, envelope modes 
|m„(^; lo)) € 7i_L, where n denotes some set of integral transverse-spatial and polarization modal indices, 
such that each mode envelope satisfies 



'i-W \'^n{£.;uj)) = H{lj) |m„(^;w)) , 



(94) 
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or equivalently 

\u^{e,co)) = u{te;^)\un{e;^)) (95) 



together with the gauge constraint H58a|l . for aU longitudinal positions —oo < £^ < oo. Because free-space 

propagation is unitary, these modes may be chosen to be orthonormal in each transverse plane, i.e., 

(u„(^;tj)|M„'(^;cj)) = (5„„/ (96) 

and complete in each transverse plane, in the sense that 

span[|«„(C;w))] (97) 

or equivalently 

"|m„(^;cc>))(m„(C;^)| =/±, (98) 



El 



where I± is the identity operator restricted to the solenoidal sub-space Ti.±, or equivalently the Hermi- 
tian projection from H into H^. Familiar and convenient choices are the usual free-space Gauss-Hermite 
modes or Gauss-Laguerre modes in paraxial optics. The spatial profiles of these modes may be charac- 
terized by specifying the longitudinal location of the focal plane and the eigenvalues of certain Hermitian 
operators of which they are eigenfunctions in that plane. For example, in the focal plane, the components 
of the Gauss-Hermite modes are simultaneous number states of two harmonic oscillator Hamiltonians 
(one for each transverse Cartesian coordinate) , while the Gauss-Laguerre modes are simultaneous eigen- 
states of a radial harmonic oscillator Hamiltonian and the longitudinal component of orbital angular 
momentum js^, H^l- Modes with slightly off-axis propagation directions can be defined in terms of 
the coherent states associated with these quadratic Hamiltonians, and many other generalizations are 
possible to suit particular problems. 

Any linear combinations of the m„ modes satisfy the homogeneous paraxial wave equation, so cannot 
possibly represent the fields in the region actually containing the sources, but they can exactly reproduce 
any fields in the "downstream" vacuum region beyond the sources, and then be extrapolated upstream. 
That is, for all ^ > £,i, the actual radiation envelope \if){(^;uj)) satisfies the homogeneous wave equation, 
lies entirely in the solenoidal sub-space H±, and may be decomposed as 

|^(e;c.))= ^|«„(e;c.))(«„(^;u;)| |^(^;c.)) 

(99) 

n 

and both the modulus and argument of the complex expansion coefficients appearing in the sum have a 
simple interpretation. Using the orthonormality of the modes, the normalized power spectral density at 
a given transverse plane ^ > is given by 

^T™[a](e; u;) = oj' {^{^ u;)^^; lj)) - ^ | {u,,{t cj)\ip{t ^)) l' ; (100) 

71 

and, individually, each 

^jTbmKK^;'^) = |(w„(e;c^)l^(C;t^))l' > o (loi) 

provides the normalized power-spectral-density contribution from the nth mode, while 

0„(^; c^) - arg[(tt„(e; u;)\rP{^; lu))] (102) 

determines the slowly- varying envelope phase (i.e., phase apart from the fast — cut dependence) of 
the the nth mode. From unitarity, it follows that ^TEM['it„](^; w), 6'„(^;w) and ^J'em[ci](?; w) = 
S ^■1'em[w„](^; w) are all independent of f for any ^ > ^i. With perhaps a slight abuse of notation, 

n 

we have used ^J'sMi'^n] with a tilde to denote the scaled power (spectral density) contribution from 
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the nth normaHzed mode i.e., the power spectral density associated through Poynting's theorem with 
the ket (m„(^; a;)|i/'(f ; w)) |it„(^;a;)) , and not the power spectral density ^^EMlMn] literally associated 
with the ket |m„(^;w)) , which is fixed at the value by our normalization conventions. 

In general, a countably infinite number of orthogonal modes are required to exactly describe the radiation 
fields, but if the modes are chosen appropriately, a finite basis set Bn oi size N, 1 < N < oo, such that 
attention is confined to n G Bn, may provide a sufficiently accurate representation. For example, if the 
mode shape is approximately Gaussian in cross section, then it should be accurately represented by the 
fundamental and perhaps a few higher-order Gauss-Hermite modes, with proper choice of the waist size 
and location. Formally, the expansion coefficients may be determined by again using the Green function 
solution (|76|) . taking advantage of the finite support of \ f{u!)) ■ In a derivation exactly analogous to that 
in the previous section, we find that for any ^ > ^i, 

min(?,?l) 



i I («„(e»i/(e';^)) ^^^^^ 



I J di'j d'r Ur.{r;^';cjr ■f{r;^';u;). 



which apart from an overall constant is just a three-dimensional inner product, or overlap integral, 
between the solenoidal source envelope f{r; and the free-space mode envelope tt„(r; w). Equiv- 
alently, we can write this as the overlap integral between the full current density j(r;^';w), and the 
normahzed electric field £n±{£,'',^) — i(jJUn{r;^';uj) associated with the nth mode: 



= -2^/ < {enA^';^)m';o^)). 



(104) 



So, again, we do not need to explicitly decompose the current density into its solenoidal and irrotational 
components in order to calculate these expansion coefficents. The factor of ^ is analogous to that 
appearing in (|87|) , and in fact H1U4|) reduces to the full energy conservation result if we are lucky enough 
to choose a basis which contains the extrapolated solution %(r;^;w) in its span. More generally, it is 
clear that the finite-basis set approximation is just the orthogonal projection of the actual extrapolated 
free-space solution into the subspace span [ |m„(^; w)) ] spanned by the modes in Bn '■ 

u;; Bn)) = ^ \u.n {(; cj) ) ( M„(e; ^)| |x(^; ^)) • (105) 



From linearity, it immediately follows that if |i/m(^;w)) is the basis-set approximation to \Xm{£,', ^)) 
resulting form a source |/m(^;<^)) then J2'^rn{(^) \i^m{£.i^)) is the approximation to J2'^rn{(^) \Xrn{S,;i^)) 

m rn 

resulting from the source envelope ^c„i(w) |/m('C;^)) , where the c„i(w) G C are arbitrary coefhcients. 

m 

So we may also determine the basis set approximation for any source by first finding the basis set 
approximation for the impulse source |r; s) and then performing a convolution over the actual source. 
That is, we can use as an approximate propagator, the restriction of the full Green operator to the 
subspace span [|«„(^;a;))] : 

|iv(e;^;6jv)) I de \un{C;^)){un{e;^)\\f{C;^)), (106) 

which clearly agrees with previous results after exchanging the order of the integral and sum. 
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It follows from the properties of orthogonal projections in Hilbert space that the projected basis-set 
approximation is also uniquely determined by the following constrained minimum-distance criteria: 



HC^^-^Bn)) = argmin || H^;lu)) - \xi^;L^)) \\ 



s.t. \v>i^;uj)) G span |M„(^;tj)) 



(107a) 
(107b) 



The estimated radiated power (spectral density) is then always a lower bound on the actual power 
(spectral density) in any plane ^ > : 

£7'^^[u]{tco;BN) = J2 £^'EMM(e;^) < £7^^^^;) = ^y.MiaKC;^), (108) 

with equality if and only if v{r; Bn) — xl'"!?!'^)- In fact, it is straightforward to establish that, 
at the constrained optimum (|107|l . the squared-distance between the basis-set approximation and the 
actual extrapolated fields is just proportional to their difference in power spectral density: 



(109) 



F. Variational Approximation 



Such connections revealed in the basis-set approach between the orthogonal projection, minimum dis- 
tance, maximum radiated power, minimal negative mechanical work, and maximal source/field overlap 
suggest that an approximate field profile for the radiation and a lower bound on the radiated power 
may be obtained directly through a general variational principle. Consider a parameterized family of 
trial envelope modes v{r; ^; lo; a) depending on some set of adjustable parameters denoted by the vector 
a, where, for any fixed choice of allowed parameter values, we assume the the trial mode satisfies the 
free-space wave equation, is solenoidal, and is normalized in any transverse plane. In principle, the trail 
mode can be written formally as a linear combination of the above basis modes: 

\vit to; a)) = J2 ^)) (110) 

n 

for some complex expansion coefficients Cn — Cn(^; ck) G C which satisfy the constraint 

^|c„(c^;a)|' = l, (111) 

n 

so that 

{viC,to;a)\vi^;Lo;a)) ^1, (112) 

but are otherwise arbitrary. Note that the expansion coefficients are here assumed independent of ^, but 
may in general depend on the adjustable parameters in some complicated, and perhaps nonlinear fashion. 
(More generally, we could allow dependence on ^ in the c„, leading to a variational problem involving 
the solution of a set of coupled ODEs rather than miniization of a function. This extra generalization is 
not really needed for the case of spontaneous radiation.) Using the Cauchy-Schwarz inequality, we have, 
for any 

I {v{^; u;; a)^^ co)) \' < {v{^; u;; a)\v{t c.; a)) o;)|^(^; u;)) , (113) 
which using (|112|l simplifies to 

(^2|(t;(e;c^;a)|V^(^;a;))|'<^2^(V'(?;c^)|t^(^;c.))), (114) 

or equivalently 

^T^MbKC;^;") < ijy^M[a]{tu;), (115) 
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with strict equality if and only if v{r;£^;uj;a) exactly mirrors the shape and polarization of the actual 
paraxial solution xp{r;^;uj) in the plane ^, up to some overall phase; i.e., 

for some real angle a). 

We have indeed arrived at a variational principle for the paraxial radiation, where we may approximate 
both the relative spatial/polarization profile and overall amplitude of the radiation fields by maximizing 
the power-spectral density in a normalized trial mode v{^;uj;(y.) measured in some post-source plane 
^ > ^1, as a function of the variational parameters a determining the trial mode's shape and polarization. 
That is, we take as the approximate radiation envelope 

Lu; a)) ^ {{v{^; lo; a)|tA(^; u;))) \v{^- to; a)) , (117) 

where the optimal parameter vector a = d[j]{Lu) is chosen so as to maximize 

^^^^(e;^^;^)^ £5',„[t;](e;c^;a) =cc>2|(t;(e;cc>;a)|t/,(e;c^))|'. (118) 

The optimized mode shape is then the best guess, within the manifold of possibilities allowed by the 
shapes parameterized by a, of the actual paraxial field profile beyond the sources, and the squared-norm 
yields a lower bound on the actual paraxial power spectral density of the radiation at the frequency under 
consideration. 

However, as written, this variational principle appears vacuous at best, since if we knew the actual 
'0(r; ^; w) so as to be able to calculate its overlap with the trial field v{r; ^; w; a), there would of course 
be no need for a variational approximation in the first place. But from H1U4|) . we may also write 

(-i;(C;t^;a)|'0(^;^)) = -2^ J d^J d^r £„i(r;^;cc;; a)* • j(r;C; w) 

5i (119) 

where e„_L(r; ^; a;; a) = iujv{r;^\Lo\a.)e'^^^ is the normalized solcnoidal frequency-domain electric field 
associated with the unit-norm trial vector potential envelope v{r; ^; lo\ a). In terms of the trial envelope 
i/(r; ^; oj; a), we see that 

£ ^.(e; uj; a) = u:^ \{v{t, lo; a)^^- a)) |' = u:^ a)\v{t, cu; a)) 

(120) 



§ / < (£.i(^';c^;a)|j(C';a;)), 



where j.(t"; w; a) = iuj i/(^; lo; a) is the solcnoidal electric field associated with the unnormalized trial 
envelope . Because the left-hand side is purely real and non-negative, the right hand side must already 
be so as well, so in fact: 



(121) 



From the inequality pi5|l and the energy-conservation result H87|l we arrive at an alternative formulation 
of the variational principle: 

£j'.(e;c.;a) < £j'^(^;^;a) ^-i^y^Je^^SjKc^ja) 

< £3'„(^;c.) = -l-y^Je^^; j](^) (122) 

= -i£j'_j£x±;j](^)- 



In addition, since 

(«.(e; cj; a)|iv(e; u;; a)) = | {v{t a)\x{^; lo)) |' = {u{i; cu; a)\x{t ^)) , 



(123) 
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for all ^, it follows that the unnormalized trial envelope i>{r;^;uj;a) (i.e., with the absolute amplitude 
included as an adjustable parameter) also minimizes the Hilbert-space distance between the trial solution 
and the actual extrapolated paraxial solution in any transverse plane ^ : 

II KC;c^;«)> - lx(e;^)>ll > II l'-(e;^;a)) - \x{^;^)) II 

The results (|122|) and H124|) are obviously closely related to the results (|108|l and (|107|l derived in the 
basis-set expansion, which first led us to the variational formulation. Equivalently, we can express our 
variational principle most succinctly in terms of the unnormalized (amplitude-included) trial wavefield 
which solves a constrained power-maximization problem for all parameters a defining the trial-solution 
i/(r;^;w;a): 

a = argmax [^T^(^; w; a)] (125a) 
s.t. £y^{^;u;;a) = -^£T^^^,[e^r,j]{io;a). (125b) 

These several equivalent formulations of the maximum-power variational principle (MP VP) are intu- 
itively appealing, perhaps so much so that they might be guessed immediately, without the entire 
build-up of mathematical machinery. The variational approximation, which maximizes the resemblance 
(in the Hilbert-space metric), within the considered family of possibilities, to the actual paraxial field 
in the post-source region, or equivalently, to the extrapolated source-free paraxial field everywhere in 
space, can be obtained by maximizing the source/field spatial overlap |/ d^' (£i,j_(^'; w; a)|ji(^';a;))| , 
i.e., the magnitude of the three-dimensional inner product between the normalized source- free trial field 
and the actual current density. That is to say, our best approximation to the free-space radiation fields 
beyond the sources is found by maximizing the physical resemblance of these fields, when extrapolated 
backward into the region of the sources assuming free-space propagation, to the profile of the sources. 

Equivalently, we can say the optimal profile is that which, if it actually were to be incident on the sources, 
would experience the maximal small-signal gain (i.e., neglecting saturation effects or indeed any back- 
action on the sources), due to energy absorbed from those sources, and in this case the virtual "gross" 
power gain (i.e., due to stimulated emission, but neglecting the stimulated absorption) is numerically 
proportional to the estimated spectral density of spontaneously-radiated power. This approximation also 
optimizes the total radiated power spectral density in the full variational radiation envelope v{r; ^; uj; a), 
consistent with the constraint, arising from energy-conservation, that this power could have arisen from 
work done by the sources, with appropriate compensation (i.e., scaling by ^) for the fact that we have 
replaced the inhomogeneous fields with their homogeneous extrapolation in the region of these sources. 

Note that this variational principle is reminiscent of the Rayleigh-Ritz variational principle familiar from 
ordinary quantum mechanics, but despite the analogies between paraxial radiation propagation and non- 
relativistic quantum mechanics, these variational principles are actually distinct. In quantum mechan- 
ics, "the" Rayleigh-Ritz method is applicable to approximating the low-lying energy eigenvalues and 
corresponding stationary states of the homogeneous, time-independent Schrodinger equation, whereas 
the present variational principle is associated with solutions to the inhomogeneous, time-dependent 
Schrodinger equation. 

However, the MPVP does share with the quantum-mechanical Raleigh-Ritz technique many of the same 
features and limitations common to variational approximations and optimization problems. Because we 
are determining a stationary point (in fact, a maximum) of a power (spectral density) functional, the 
actual value obtained for the power spectrum at the extremum is relatively insensitive to the precise 
shape of the trial field ~ generally, if the characteristic relative error in the latter is 0{e) in some sense, 
then the relative error in the former is O(e^), because the first-order variations vanish. In order to 
estimate the power to say 1%, we need only match the variational parameters describing the field shape 
to about 10%. This is, of course, great news if we are actually most interested in estimating the power, 
but, conversely, much less satisfying if we are interested in determining, say, the spot size wq, or some 
characterization of the the angular spectrum. This is analogous to the situation in quantum mechanics, 
where it is well known that the Raleigh-Ritz approximation estimates the ground-state energy eigenvalue 
more accurately than its corresponding wavefunction. 

The variational approximation for the field profile is expected to be accurate to the extent that the 
trial mode can overlap, or mimic, the actual mode in any transverse plane ^ > in the free-space 
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region beyond the sources, and obviously the approximation is expected to improve as the number of 
functionaUy-indcpendent shape parameters is increased. Since lu; a) satisfies the free-space equation, 
it is actually an approximation throughout all ^ to the extrapolated envelope x{''"t£,!<^)- The correspond- 
ing power yields a lower bound for the actual power radiated, but in the absence of an upper bound, one 
does not have any foolproof means to estimate the closeness of this approximated to actual power. We 
can consider a so-callcid "minimizing" sequence of variational problems where at each stage we expand 
the parametric family of trial envelopes to include additional possible details of the field profiles which 
we anticipate might be present in the actual radiation but which could not be captured in the previous 
trial fields. If the extended family continues to include as a possibility the previous variational solution, 
then the estimated power monotonically approaches the actual (unknown) power from below, and the 
Hilbert-space distance between the trial envelope and the actual (unknown) envelope in any post-source 
plane monotonically approaches zero from above. In practice, we can stop when the difference (in power 
or Hilbert-space distance) upon adding finer details to the trial envelope becomes sufficiently small, or 
else when we have expended as much computational effort as we can spare. From above, we know that 
an additional parameter leading to an 0(e) relative increase in the power is expected to produce only 
an O(-ye) relative improvement in the local field profile. 

Although we formally expressed the normalized ket \v(i^;uj;a)) as a sum over the orthonormal modes 
|m„(^; (jj; a)) , we stress that the variational envelope can be c;xplicitly written in any arbitrary form, with 
adjustable parameters which appear linearly or non-linearly, as long as it satisfies the homogeneous wave 
equation, gauge constraint, and normalization constraint. If the parameters a represent the expansion 
coefficients in some finite linear sum of orthonormal modes, which is the assumption of the c;lassic 
Raleigh-Ritz and Ritz-Galerkin techniques in the calculus of variations, then the variational approach 
formally reduces to the basis-set approach, but it is often more convenient or efficient in practice to use 
an explicitly nonlinear family of trial solutions where the variational parameters do not appear linearly. 
For example, one might use a Gaussian beam with indeterminate waist size and location, as discussed in 
the introduction, or perhaps even additional parameters to include skew or eccentricity or higher-order 
structure in the spot profile. 



IV. NON-PARAXIAL GENERALIZATION 

When translated from the quantum mechanical into conventional notation, none of the energy-balance 
and variational results derived above in the paraxial limit depend in any essential way on the assumed 
paraxial nature of the fields. This suggests that analogous results should hold more generally, beyond 
the paraxial approximation, for any radiation emitted spontaneously by charges moving along prescribed 
trajectories. The essential features which emerged in the paraxial case, which we now abstract to more 
general problems, are: approximation of the actual radiation fields arising from the sources by homo- 
geneous (i.e., free-space, or source- free) fields; a variational principle, derived from the Cauchy-Schwarz 
inequality in a Hilbert space picture, mandating the maximization of (only) outgoing power spectral den- 
sity: c;nergy exchange between fields and charges which can be characterized as the three-dimensional 
inner product, or overlap integral, between the extrapolated, solenoidal electric field and either the full or 
solenoidal current density, as convenient; and, finally, a functional relationship between the mechanical 
work performed on/by the charges by/on the extrapolated fields within a region and the radiation flux 
through a boundary surface, arising, physically, from energy conservation, and, mathematically, either 
explicitly from the formal solution or else from some form of the the Fundamental Theorem of Calculus 
(e.g.. Gauss's theorem or Green's Identities in multiple dimensions.) 



A. Vector Spherical Harmonics and the Spherical Wave Basis 

In order to generalize our results to the non-paraxial case, we first introduce a bit more mathematical 
notation and machinery. We use the scaled spatial position ^, and introduce spherical coordinates: the 
scaled radius C = ||CII ^ 0) the polar angle 6, and the azimuthal angle ip, so that C = CC = C^- We let R 
denote any (closed) spatial region in R^, and dR will denote its boundary. For any position and any 
nonnegative (' > 0, we define the closed ball of radius (' centered at Co as i?(C'; Co) = {C I IIC ~ Co|| < C'} > 
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and we will take V(C';Co) C IR^ to denote a closed, simply-connected region satisfying B{(';Cq) C 
^o)- For simplicity, whenever the central position coincides with the origin {Cq = 0), the dependence 
will be suppressed in the notation; i.e., B{C) = -B(C';0) and V{C) = ^^(C';0)- 

As in the paraxial case, we take the (scaled, frequency-domain) current density w), assumed known, is 
to be localized in space, vanishing identically outside some finite-radius ball B{([), for some < < oo. 
Again, the solenoidal component J_l(C;^) will therefore not have have compact support in general, but 
because of the rapid falloff of its spatial dependence (0(C^^) as C oo), we may, with arbitrary small 
error in the calculated radiation fields, neglect it beyond some sufficiently large but finite radius Ci ^ Ci- 
(Again, at the end of the calculation, we could attempt to take this support to be infinite, provided that 
the full current density continues to fall off sufficiently rapidly to ensure convergence of the integrals, 
and to enable somehow the unambiguous characterization of the "far-zone" versus non-radiative fields, 
when the latter extend to infinity.) 

It will prove very convenient to have some explicit basis in which to decompose the vector potential. 
The plane-wave (spatial-Fourier) basis immediately comes to mind, but certain cumbersome singularities 
inevitably arise in this representation. Namely, whenever the vacuum dispersion relation is satisfied, i.e., 
1 1 A; 1 1 = Lo, solutions to the homogeneous Helmholtz equation possess delta- function singularities, while 
solutions to the inhomogcneous solution exhibit second-order divergences whenever the corresponding 
Fourier coniponents j±{k,Lu) ^ 0. In order to avoid such difficulties, we will instead use the spherical 
wave basis 0, HE ■ The conventional vector spherical harmonics are defined as 

Xg,n = ^£m(C) = Xi„i{9, if) = --j=k=LYem{d, f), (126) 

where, in analogy to quantum mechanics, we have defined the Hermitian (scaled) orbital angular mo- 
mentum operator as 

L = (53D X = C X = (C X a) , (127) 

which acts only on the angular degrees of freedom (i.e., not on the radial or the polarization degrees-of- 
freedom) , 

^ = 72^+^'^" (-M +*cot(0)^ ) + -^^e^e+^^ (+^ +^cot(0)^ ) - '^W' 

and the functions Ypm{9^Lp) denote the usual scalar spherical harmonics for nonnegative integers t — 
0,1,2,..., and integer m = — £, — ^ + 1, . . . , 0, . . . , ^ — 1, ^. In addition, we define two other related sets 
of basis vector fields: 

Ze.m — ZimiC) — Zgmid,'^) = C X Xg^, (129) 

and 

Nim = Ni^{C)^Ni^{e,ip)=CYem. (130) 

(By definition, we take Xqq — Zqq = identically, ultimately reflecting the fact that spherically- 
symmetric solutions to the free-space Maxwell equations can exist only in the non-radiative a; — > 
limit.) 

The orbital angular momentum operator L effects infinitesimal (inverse) rotations of the spatial degrees 
of freedom, i.e., observation position or field point, and satisfies the useful identities: 



C-L =0, (131a) 

LxL^iL, (131b) 

[L^,L] = 0, (131c) 

[Pl,L]^0, (131d) 

L' = r-^{r )+r^Pl. (131e) 



The scalar spherical harmonics themselves may be interpreted as simultaneous eigenstates of and Lz ■ 
Each of the basis vector fields may then be interpreted as a simultaneous eigenstate of = + Jy + 
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and Jz for some spin-one object, where J = L + S \s the total (spatial plus spin/polarization) angular 
momentum operator. For fixed I and fixed W = X , Z , or iV, the members Wirn{^^ f) of each of these 
three families therefore transform amongst themselves under total spatial rotations (of both spatial and 
polarization degrees of freedom). The spin angular momentum operator S generates rotations of the 
polarization vector, and in the Cartesian basis, S = xSi + yS2 + zS^ may be defined as a vector of 3 x 3 
Hermitian matrices Sj for j = 1, 2, 3, whose components are given by 

where Cjkn is the completely antisymmetric Levi-Civita tensor in three dimensions. Prom this definition 
one can prove the useful operator identity 

dx =P3D-S. (133) 



Using various vector identities, these spherical vector harmonic basis fields also can be shown to satisfy: 



and 



as well as 



and 



and also 



together with 



L-Xem^ V£{iTT}Yem, (134a) 

L-Zem = 0, (134b) 

L-Ne^ = 0; (134c) 

L xX^rn = iXe„i, (135a) 



c 


■Xirii 


= 0, 


c 


■ 


= 0, 


c- 




= Ye 



LxZim = VW+^Nem, (135b) 
LxNem= 0; (135c) 



(136a) 
(136b) 
(136c) 



CxXirn= Zirn, (137a) 

C X Zf rn = ~^{m, (137b) 

CxNe^= 0; (137c) 

0, (138a) 



d 


■Xim 


d 


■ 


d- 





= -^VW+^)yem, (138b) 
= ^Yem-, (138c) 



dxXem= i {i^iie + l)Nem + Zem) , (139a) 

dx Zi^^-^Xim, (139b) 

dx Ne^^-^^eie+l)Xi„,. (139c) 

For any allowed £ and m, the vector spherical harmonics are geometrically orthogonal as vectors at every 
point ^ on the unit sphere: 

XemiCr ■ ZtmiC) = Xt^iO* ■ NemiC) = ZtM* ■ NemiC) = 0; (140) 
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are (Hilbert-space) orthonormal when integrated over solid angle (except in the trivial case where both 
vector fields identically vanish for £ = 0): 

j dO sin 6* J dip Wi ,n{0 , ip)* ■ W^,„^,{9,(p) = 5ww'Sii'S,n„i' (1 - SwxSio) (1 - SwzSeo) ■ (141) 
and also collectively constitute a complete set for vector fields defined on the unit sphere: 

oo I 
W=X,Z,N £=0 m=-e 

where I3 is the identity matrix on the polarization degrees of freedom. 

Now, outside the (effective) support of the solenoidal sources, i.e., for any ( > (,2 > Cii the actual 
Coulomb-gauge vector potential a(C;w) is assumed to satisfy the homogeneous Helmholtz equation, 
and the method of separation-of-variables may be used in spherical coordinates to show that there the 
solution may always be written in the form: 

a(C;c.) = J2 ^^^Md X [/i«(fcC)^£™(^?, V')] - t<^r^iioW^\kOX,„,{9,^), (143) 

I, m 

where the spherical Hankel functions of the first and second kind, 

h''^\x) ^ ji{x) +ineXx), (144a) 
hf\x)=Mx)-ineix). (144b) 
are defined in terms of the usual spherical Bessel functions, 

1 d \ ^ / sin X 



^^^^^ = ^-^y[-xd-x) l^j' (14^) 
(which should not be confused with any current density), and the spherical Neumann functions, 

/ ^ d \ ^ / cos a; \ , , 

ne{x) = -i~xYi-—) . (146) 

\x dx / \ x / 

Any two of these four families of special functions may be taken as the linearly-independent solutions 
to the radial differential equation obtained through separation of variables. The appearance of only the 
h^^\kr) (X in the radial part of our formal solution reflects the facts that the solution holds in the 
extra-source region C > Ci which excludes the origin, and that we have imposed the so-called Sommerfeld 
boundary conditions, requiring only outgoing radiation at spatial infinity. 

Each term in the superposition H143() is everywhere divergenceless, so the total vector field is solenoidal, 
as required. The expansion coefficients a|^(a;) G C and a"^(a;) € C are, respectively, the (scaled) 
electric and magnetic multipole moments of the transverse current density j_l(C;w), and are uniquely 
determined by quadratures over j, or equivalently by field values in the far-field (or, in fact, by just the 
radial components of the fields on any two spherical shells outside the sources.) For example, using the 
far-field asymptotic behavior, one can, with a bit of algebra, show that 

max[m+l,. 



Q"mH = - ^ E M,,r^'{Lo)- jdes{nejdp,Yi„,{e,^)*LYt„AO.^), (147a) 

V ^ ~^ m'— min[0,m— 1] 
max[m+l,£] 

4^{u)= *^|^=y E Rira'{Lo)- J d9 sinej d^Ye^{d,^rLY„r,,{9,^y, (147b) 



where we have defined 



m'— min[0,m— 1] 

Mernicj) = J d^C [Mkr)Ye,nify]j±{C;i^), (148a) 
= / d^C [jiikr)Y,Urr]dxj^{C;Lu). (148b) 



26 



Many other equivalent forms are possible, depending on precisely what data are used and where, but 
explicit forms for the multipole expansion coefficients will not actually be needed here. (After all, if 
they were known or easily calculable, we would not need to resort to any sort of approximation.) Note 
that time-reversed, or equivalently conjugate, or incoming- wave, solutions can obviously be written in 
an analogous manner, with the (kC) replaced with (^C) c< ^tq~- 

By assumption, any solenoidal, source-free "trial" vector potential x{C'i^'iOl) satisfies the homogeneous 
equation everywhere in space, and may be expressed, for any C,, in a similar form: 

X(C; ^; a) - ^ ^Xl„,{^; a)9x [j,(fcC)X,„(0, ^)] + -^x7,^{oj- cx)jdkC)Xe„,{0, ^). (149) 

i, rn 

The use of je{k() is the only choice for the radial dependence which solves the Helmholtz equation but 
is regular everywhere in space, including the origin. Each term in p49|l is also automatically diver- 
genceless everywhere, so their sum (when convergent) is solenoidal. The complex coefficients Ximi'^^ ^) 
and xtmi^:^) for this source-free solution may nevertheless be interpreted as some effective (scaled) 
multipole moments which would describe what we might term an "effective source" g{(^;oj; a), or may 
just be regarded as free expansion coefficients determined by (or even taken to be equal to) the actual 
variational parameters a appearing in x(^;a;;Q:). The familiar expansions 

oo £ oo 

= AnY^^'MkC) J2 ^'^™(C)*>£m(fc) = J] i'VM'^^ + l)j£(fcC)Fi o(arccos(C • fe),0), (150) 

1=0 m=-e e=o 

for any fixed (scaled) wavevector k — kk, may be used to establish the connection between the plane wave 
and spherical wave representations, and to confirm that they both span the same space of source-free 
solutions. 



B. Green Functions 



Alternatively, the exact solution to the inhomogencous problem at any position ^ may be expressed in 
terms of a convolution over the retarded Green function: 

+4fe||C-c'|| 

which is invariant under simultaneous translations, rotations, or reflections of both spatial arguments, 
and satisfies the scalar Helmholtz equation with impulsive source (and with negative unit charge, con- 
forming to convention for electromagnetic problems), 

(92 + c.2)G'-(C;C';^) = -'5(C-C'), (152) 

together with the Sommerfeld asymptotic radiation conditions for outgoing waves: 

G"'(C; C; ^) = 0(1/0 as C -> oo (153a) 
(^-zA:)G'-(C;C';^) = o(l/C) as C^oo, (153b) 

for any fixed emitter position (Because we are working with solenoidal sources in otherwise free space, 
we can use a scalar rather than dyadic Green function [sj.) The retarded Green function G"'(C; C; 1^)1 
then represents an outgoing (scalar) spherical wave at the observation position ^ emanating from a 
harmonic point source at (By symmetry, analogous conditions and interpretations obviously hold 
true under permutation of the two spatial positions C a-nd C^ from which follow various well-known 
electromagnetic reciprocity results such as the Raleigh-Carson and Lorentz theorems.) The advanced 
Green function 

G-^(C;C';^) = [G-(C;C';^)] = a^c-CT ^^^^^ 
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satisfies the same impulsive Helmholtz cqiiation, and in fact may be interpreted as the time-reversal of 
the retarded impulse response, representing an incoming spherical wave, observed at ^, and ultimately 
converging to a point absorber at Given any current density j{C';i^) with compact support in the 
neighborhood V{(i) of the origin, and supposing no incoming radiation from distant sources by assump- 
tion, the solenoidal component J_l(C')'^) can result only in outgoing radiation asymptotically as ( ^ oo, 
so the Coulomb-gauge vector potential can always be expressed as: 

a{C;io)= j d'CG-\C;C;^)j±{C;co)= j d'C G-'{C;C;u>)j±{C;u>). (155) 

viCi) 

Clearly, this satisfies the appropriate inhomogeneous Helmholtz equation: 

{d^ + oj^)a{C,u;)= J d^' {d^+u;^)G"^{C,C';uj)j±{C';uj) 

'5(C - C) jl(C», = (C; cc) 



(156) 



and also shares with G'°*{<^; C;uj) the appropriate outgoing- wave Sommerfeld boundary conditions. Be- 
cause G""(C; C';'^) is radially symmetric, depending on the spatial coordinates only through \\<^ — ^'||, 
and because d ■ j± w) = by construction, we have 

d- [G-(C;C^'^)i±(C^'^)] =SG-(C;C^c^)-j±(C^c^) = -a'G-'(C;C^'^)-i±(C';'^) 

= -d'-[G-^{C;C;u^)3±iC;u^)]- ^ ' 

Since G"*(C;C';w) = 0(1/ |C - C'|) and j_L(C;w) < 0(1/C^) as C ^ oo, Gauss's law then implies 

d-a{C,u^) = j d\'d-[G-^{C,C,';o^)3±{C,';^)\ = - j d\'d'-[G-\C-X'\oj)3^{C'\oj)] 

f , , (158) 

= - hm / da'h ■ [G"*(C; C'; ^)3±{C; <^)]=^ 

C2^00 J 

so indeed a(^;w) is solenoidal, as required. (Subsequently, we will make use of a number of similar 
multidimensional "integrations by parts," relying on the the sufficiently rapid fall-off of j± and radiative 
fall-off of G"*, but for brevity we will not always provide the explicit demonstrations.) 

Outside the support of sources, say for C ^ C2 > Cij the solution to the inhomogeneous Helmholtz 
equation may also be written in terms of a surface integral involving the Green function and boundary 
data. Specifically, by applying Green's second identity to the region lim ^^(Cs) — V{(2), using the 

Helmholtz equations satisfied by the vector potential and the retarded Green function, and making use 
of the assumed Sommerfeld radiation conditions, we find: 

a{C,iv)= J da'[a{C;u){h'-&)G-\CX'-,^)-G-'{CX'-,co){h'-d')a{C-,u;)], (159) 

dV(C2) 

which is the vector form of the well-known Kirchoff diffraction integral. Here the unit vector n' = n(C') 
denotes an outward normal to the boundary surface at any point C G dV{(2)- 

By piecing together different spherical wave solutions to produce just the right discontinuity at C = C) 
it is straightforward to establish the well-known Bessel decomposition of the retarded Green function: 

00 t 

G-^{C,C';oj) = iojY,MkC<)h'-l\kC>) YemiOYemiCr (160) 
e=o m=-e 

where 

C< =inin[C,C'], (161a) 
C> =max[C,C']. (161b) 
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From linearity, we can see that the difference between the retarded and advanced Green functions, 

DiC; C';^) ^ G-(C; C';^) - G--(C; C; ^) - = ^ ^^^^ (^IK - C'll) , (162) 

must be a solution (with respect to either spatial coordinate) to the source- free scalar Helmholtz equation 
everywhere in space. Substituting in the decomposition IjltiOf) for the Green function, this homogeneous 
solution D{C; C'','^) may be written as: 

D{C;C';u;) = 2iu;J2MkOMkCrYi„,{C)Ye^{Cr, (163) 

^ . . 

where we have used the facts that ^ Y(m{C)Yem{C')* f^nd je{kC) are both real. 

Now, it turns out that, just as a general solution to the inhomogeneous equation with outgoing Som- 
merfeld radiation boundary conditions may be written as a convolution between the retarded Green 
function and the solenoidal portion of the current density, any arbitrary solenoidal solution to the homo- 
geneous problem may be written as a convolution between what we may term the fundamental solution 
D{(^; oj) and the solenoidal component g±{C'', t^), of some effective source term g{C', ^) whose support 
can be confined to any finite-radius ball i?(Co) where > : 

x(C; ^) = y d'C D{C; C; ^)s±(C'; = / ^'C d{C, C; ^)9±(C'; ^) (164) 

V{Co) 

To see this, first note that if we choose 

giCiu) = g^iCiu) ^ e{Co ^ OjtikOXirniO,^), (165) 

(which is clearly divergenceless by the chain rule and the transverse nature of XirniO)^ then, by using the 
explicit decomposition (|163|l for D{(^; uj) and the orthogonality properties of the spherical harmonics 
as well as their behavior under the action of the orbital angular momentum operator i, we find 



Co 



MkOXemie^f), (166) 



and thus any term of the magnetic multipole form in (|149|l can be generated in this manner. Using 
standard vector identities, the symmetry of D(^; lu) in its spatial arguments, and integration by parts 
relying on the rapid fall-off of g± , we also have 



dx J d3^'i?(C;C';^)9(C';^)- J d^c'ax [z?(C;C';^)s±(C';'^)] 

= J d'CD{C,C;^)d'xg^{C;u;), 



(167) 



So any required term of the electric multipole form can also be so generated. By linearity, an arbitrary 
homogeneous solution may then be constructed as a convolution over the fundamental homogeneous 
solution. We will see below that the effective source g{C;uj) appearing in the convolutional representa- 
tion may be taken to be identical to that determining the multipole coefficients in the spherical wave 
expansion. 



C. Hilbert Space Results 



In the non-paraxial case, the wave-equation is no longer equivalent to a time-dependent Schrodinger 
equation, although the source- free Helmholtz equation is analogous to a time-independent, free-particle 
Schrodinger equation, albeit typically with non-normalizable, scattering-state boundary conditions. 
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However, when convenient, we can still recast our non-paraxial mathematical results into a Hilbert- 
space-like formalism. Without a distinguished optic axis, we extend the spatial dependence of spinors 
and observables from two dimensions to three dimensions, each coordinate now treated on an equal 
footing, and we allow for more general polarization bases, using some given family (labeled here by n, 
standing for some set of integer indices or "quantum numbers") of smooth, possibly complex- valued 
vector fields e" (C;w) (still indexed, say by s = — 1,0, +1) which constitute a local orthonormal triad at 
every spatial position (J : 

e:(C;a;)*-e^,(C;a;)=<5,,,. (168) 

For example, we could use, as convenient, a fixed Cartesian basis x, y, z: the spherical-coordinate basis 
6, (p; the helicity basis e+, e_, z (also known as the spherical basis because of its role in defining 
spherical tensors, not to be confused with the above basis vectors for spherical coordinates); normalized 
vector spherical harmonic basis vector field for some £>1 and m, i.e., Xi„i{C), N'^m(C)j ^em{C)] or a 
natural radiation basis field derived from some reference vector potential a{(^;oj), namely ia±, •So, ba- 

Much of the "quantum" -like formalism developed in the paraxial case can be generalized in the obvious 
way to a three-dimensional Hilbert space Ti^o consisting of square- integrable functions g : M.^ C^, 
and will not be reproduced here, except to note that we will now denote the (possibly generalized, i.e., 
non-normalizable) eigenket corresponding to the vector field g{<^; cu) by \g{co)) , and the generalized 3D 
position/polarization eigenkets by |^; e"(^;a;)) , such that 

{C;e^{C;oj)\g{uj)) = e?(C;c.)* • g(C;a;), (169) 

where the Hilbert-space inner product now involves integration over all spatial variables; 

{9W)=Jd\g{Cr-g'{C). (170) 

Now, consider the Hilbert-space operator K = K{lu) defined via its spatial/polarization kernel: 

K{C, s; C, s'; w) = (C; e:(C; u,) \K{u)\ C; e^, (C; a;)) = -^iuj D{C, C», (171) 
which is given explicitly by 

i^(C,0;C^0;a;) = ^sinc(fc||C-C'll)=c^'EJ■^(^Oi^(^C')*^^™(^)^^-(^')*• (l^^) 

e,m 

The kernel is seen to be real- valued and symmetric in all spin/spatial arguments, so the corresponding 
operator K ^ K'^ is Hermitian. Because the kernel acts as the identity on polarization degrees of freedom, 
the operator is independent of spin observables, and may be formally written as a function only of spatial 
observables. In fact, because the kernel is translationally and rotationally invariant, depending on the 
spatial coordinates only through ||^ — ^'||, K must be a function solely of the momentum magnitude 
II-PsdII = V—d^- Crucially, K is also positive semidefinite, in the sense that 

{g{iv) \K{co)\g'{u,)) = jd'cj d'C if(C,0;C',0;a;)g(C;w)*-9(C';w) > 0, (173) 

for any complex vector field g{(^-,uj) which is sufficiently well-behaved for the integral to exist. This can 
be seen most easily from the diagonal expansion in spherical waves: 



{g{uj) \K{u;)\g'iu;))=uj'Y.Yl I '^'^ ^:(C)*-9(C;'^)i^(fcC)*:^^m(C)' 



> 0. (174) 



It is also a direct consequence of the well-known Bochner theorem, which says that a continuous function 
in real-space is nonnegative definite as a spatial kernel if and only if its Fourier transform is proportional 
to a positive measure over reciprocal space (and positive definite if this measure is everywhere non-zero). 
Transforming to the scaled momentum (i.e., wavevector, or reciprocal-space) representation where the 
operator K is diagonal, we find: 

K{UJ) = 7r5(P3n-P3D - ^2), (175) 
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so that 

(gH|ifM|g'H)=7r J d^k |g(fe; c.)!^ 5(||fcf - c.^) > 0, (176) 

where k is the scaled wavevector conjugate to ^, and g{k\Lo) is the (scaled) spatial Fourier transform 
oi g{C,\Lo). Note that K{Ld) acts to restrict the spectral content of any vector field g{C,\Lo) in its domain 
to the manifold of wavevectors satisfying the vacuum dispersion relation lo^ — ||fc|p. However, K is not 
idempotent, so cannot be considered a true projection. 

That K{C,\C,''t^) is only positive semidcfinitc and not strictly positive definite is a consequence of the 
well-known existence of non-radiating sources for the Hclmholtz problem, which invariably complicate 
uniqueness results and inverse-scattering calculations. Specifically, suppose we take 

go(C;^) =so±(C;^) «e(Co-C).9£(fcC)^^m(^,<p), (177) 

for any allowed I and m, some > 0, and any well-behaved, non-trivial scalar function gi{kC,) which is 
chosen to satisfy 

Co 

j dCegt{kOji{kO=Q, (178a) 





Co 



dCCm{k09i{kO>^, (178b) 







Then using the spherical wave expansions, it immediately follows for this source that, for all 

xo(C; ^) - y d?C' D{C, C; ^)go(C'; cj) = o, (179) 

and also, anywhere in the extra-source region, i.e., for C > 

ar (C; cj)= J d^C G"'(C; C; cj)go{C'; ^) = O, (180a) 

<(C; c^) - / d\' G"^^(C; C'; ^)9o{C'\ u^) = o, (I80b) 

and indeed no radiation is produced by go acting either as a real, time-reversed, or effective (homoge- 
neous) source. Any such vector field go{Ci ^) corresponds to a ket \go{uj)) in the nuUspace of the operator 
K; i.e., K{u;) \go{io)) = 0. 

The nuUspace of the Hermitian operator K is therefore infinite-dimensional (and in fact non- 
denumerable). Modulo some (admittedly thorny) normalization issues due to the unbounded nature 
of K, the basis of homogeneous spherical waves solutions H149() are the orthogonal eigenfunctions of K in 
the position basis corresponding to non-zero eigenvalue. Although it would involve more mathematical 
detail than justified in this presentation, the properties and behavior of K can be analyzed rigorously 
within the theory of Reproducing Kernel Hilbert Spaces (RKHS), where the spherical wave expansion of 
the kernel emerges as a consequence of Mercer's theorem jSSi] . In a similar manner, we may define Green 
operators G'"*(tt') and G'"'"'(lj) from their spatial representations; these operators are complex-symmetric, 
but neither Hermitian nor positive semidefinite. 

However, the problems encountered with normalization cannot be completely ignored. While this Hilbert 
space TisD, equipped with the Euclidean/L2 inner product, is the natural setting for the source-terms such 
as j(Ci'^)) J-L(Ci'^)) 9{Cl^)j ^nd g±{(^;u!), it is not always a comfortable home for the vector potentials 
a{(^;ui) or xiCl'^) oi the corresponding fields, since the presence of far-zone radiation fields, falling 
off like 0(1/C), implies that the vector potentials and fields are not normalizable with respect to this 
inner product, so do not strictly belong to the Hilbert space Hso- However, they may be approximated 
arbitrarily closely (point-wise) by vector fields in the Hilbert-space, and inner products between radiation 
and sources, e.g., {j {uj)\a{Lu)) will exist for allowed current densities, which must be localized. 
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The lack of nomalizability may be traced to the fact that radiation fields can contain, in principle, an 
infinite amount of energy when integrated over all space, while transmitting only a finite amount of 
power. If necessary, we could use the familiar regularizing device of introducing a finite "quantization" 
volume V{(^q) together with periodic or conducting-wall boundary conditions, and then take the limit 
— > (X) at the very end of the calculation. For the currently envisioned application, this is not terribly 
convenient, as we are interested primarily in the far-field from the start. Alternatively, we can choose 
a different inner product, which effectively normalizes the radiation fields based on power (or power 
spectral density) rather than energy. We can define a Hilbert space Ti^^t spanned by the basis consisting 
of outgoing solenoidal vector spherical waves, with an inner product defined as a Euclidian dot product 
on the multipole expansion coefficients: 

(a(c.) , a'(^))„^^ = ['^Trni^HlA^) + aZr{^)ar„A^)] . (181) 

The exact relationship between this inner product and the scaled power spectral density will be seen 
shortly, but it should be apparent that this quantity may be expressed soley in in terms of far-field 
data, rather than on the field throughout all space, and physical radiation fields, with finite power, 
will be normalizable in the sense (a(w), ^•(w)) < oo. An isomorphic Hilbert space 7ii„ = Ti*^,^ = Tiout 
may be defined in an exactly analogous manner for incoming-wave solutions, such as a(C]Lu)*. Complex 
conjugation provides a natural isometry between them, or instead one may use the mapping which simply 
swaps the spherical Hankel functions /i*-^^ (fcr) and ft,^^^ (kr) in the spherical wave expansions without 
changing the multipole expansion coefficients or the vector spherical harmonics. Strictly outside the 
support of the sources (if any), any solenoidal solution xp to the vector Helmholtz equation (homogeneous 
or inhomogencous), of finite power but with otherwise arbitrary boundary conditions, may be expressed 
as the superposition of incoming and outgoing spherical waves, so it necessarily lies in the Hilbert space 
Tiout ® '^in- Thus, we can in this case uniquely decompose (for <^ > (2) any such solution tpiC', t^) as 

^(C;c^) = iA""(C;^) + t/''°(C;c^), (182) 

where ■0°"' and are, respectively, the outgoing and incoming components of the vector field, in the 
sense that they satisfy the appropriate Sommerfeld boundary conditions as C ^ cxd. The inherited inner 
product is given by 

(t/>(u;) , xP'iu;)) = (t/,°-(u;) , t/,'™'(c.))„,, + (t/>'»(c^) , t/,'"^(^)).„ . (183) 

Our approximating "trial" vector potentials x will lie in the proper Hilbert subspace Tihom ^out ffi^im 
consisting of source-free solutions which are regular everywhere in space. This space of homogeneous 
solutions is also isomorphic to the outgoing or incoming wave spaces: Tihom — ^out — Win, as can be seen 
by imagining the mapping that replaces, in the spherical wave expansions, each real-valued spherical 
Bessel function ji{kr) with the corresponding complex spherical Hankel function h'^^\kr) or h'''^\kr), 
respectively, without altering the angular dependence or the expansion coefficients. 



D. Energy Balance and Poynting Flux 



Using the Helmholtz equation, gauge constraint, and various vector identities, one can derive a general- 
ized impedance relation, or complexified power-balance equation, for the fields in the frequency domain 
over any arbitrary spatial region R, and valid for any consistent boundary conditions imposed on the 
radiation at spatial infinity: 

Jd^C[jAC:i^r-eaAC,Oj)]+lLoJd\ [|l£ai(C;^)ll'-|lba(C;^)ll'] + J d'<J [u- S^{C, io)] = 0, (184) 
R R dR 

where h = n(C) is a outward normal unit vector on the boundary surface at € dR. This relation is 
quite similar to the usual conservation law derived in textbooks from the time-harmonic formulation of 
Maxwell's equations 2, 5], except that it involves only the solenoidal components of the sources and the 
fields - the longitudinal component j|| of the current density and the scalar potential and its derivatives 
do not appear. 
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In the most general, inhomogeneous case, this relation implies: 

-Re J £C [jl-£a±]=Re j [fi-s^], (185) 



dR 



which is an expression of time-averaged energy conservation, relating the (negative) rate of mechanical 
work done on the solenoidal sources within R to the rate of energy escape in the form of radiation 
through the boundary dR. If we let R = V{(2) for any ^2 > Cij '^e can replace j± with the the full 
current density j on the right-hand-side, and this becomes 

-£?_Je„i;j](cc;) = -Rcy"d='C [j*-£„J=Rc J d^a [h-s^] = £'P^^,[a]iC2;uj) (186) 

dV(C2) 

The radiated power spectral density is independent of radius C2 as long as C2 > Cii as assumed, so we 
may take C2 ^ 00 on the left-hand-side whenever convenient. Using Green's Second Identity and other 
vector identities, the spectral density of the radiative flux can also be written in the following useful 
forms within the Coulomb gauge: 

^T,M[a](C2;^) = y y"d\ [a.vV-a*.V2a] =y J d^a [a-{n-d)a* - a* ■{ffd)a] . (187) 

V(C2) dV(C2) 

Similarly, if a{(^;uj) is the exact solenoidal radiation corresponding to j±{C;uj), and xiC^^) is any 
arbitrary solenoidal solution to the source-free Helmholtz equation, then for C2 > Cij 



j d-a [x-(n.a)a*-aMn.a)x + a.(n.a)x*-xMn-a)a] 



(188) 



For any such solution x(C;'-^) to the homogeneous Helmholtz equation (i.e., where d^x — ~'^'^X)^ the 
Lorentz invariant [||£x-lIP ^ ll'^xlP] vanishes everywhere (as can be seen easily, for example, from an 
expansion in transverse plane waves, or from the spherical wave basis), and from H184|l or (|188|l we have: 



1 [n-s^] = 1 d\ {d-s^\ = y/rf'C [X-5V-X*-9'X] =0. (189) 

dR R R 

That is, both real and imaginary parts of the Poynting flux vanish when integrated over any closed 
surface, reflecting the fact that, as is intuitively evident from the convolution representation (|164|l or 
from a plane-wave decomposition, any vector potential which satisfies the homogeneous wave equation 
everywhere must, in fact, result in an equal amount of averaged radiative flux passing into and out from 
any closed surface, without any dissipation of field energy or reactive storage of energy in localized, 
quasi-static fields. 

As we saw above, however, we can locally decompose the homogeneous fields into what represent, or at 
least eventually (as C ~^ 00) become, separately outwardly-radiating and inwardly-radiating waves: 

x(C;'^) = x°"'(C;w) + x'"(C;^), (190) 

where for all C G dy{C,'i)^ Re [n(C) • ■s^"'] > and Re [n(C) • s^'] < as C2 ~> 00. This is most easily ac- 
complished simply by separating the fundamental representation into retarded (outgoing) and advanced 
(incoming) parts: 



x°"'(C; ^) = y d'C G"'(C; C; ^)9±(C'; 



(191) 
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and 

x-(C;t^) = - J d'C'G-^(C;C';w)9±(C';'^)- (192) 

Note that these vector potentials separately satisfy the homogeneous equation only for C > Coi but 
satisfy the inhomogeneous equation everywhere for the appropriate effective source terms and specified 
boundary conditions. Although X'"{C]^) involves the time-reversed (advanced) Green function, it is not 
the time- reversal of X°"'(C; w), because it involves the negation, not complex conjugation, of the effective 
source g±{C', <^)- Indeed, the conjugate wave X'"(C; ^)* ^s-Y be interpreted as the usual outgoing radiation 
from the effective source —g±{(^';uj)*. 

Using the symmetry of the Green functions and a simple change of variables, it is easily shown that for 
any C2 > Co, 



d 



= (g^{Lu)\K{u;)\g^{Lu))>0, 



(193) 



which vanishes if and only if g_\_{C'j ^) is a non-radiating source at frequency lo. In the same manner, we 
also find 



(194) 



£3'bm[x"'](C2;^) = -^T_„[£;^^;-g](c.) = -J d'c J A' Kicx';^) 9±{C,u^r ■g±{C ■,u;) 

^-(g^{cu)\K{oj)\g^iLu)) <0, 
also vanishing if and only if gj_{C;Lu) is a non-radiating source. Thus, we see that, as anticipated, 

|jJ'bm[x°"'](C2;c.) = -^T™[x'1(C2;c.) - ^J'.m[x"1(C2;c.) > 0. (195) 

It then follows that 

^ijd'cj d'C zc.i?(C;C';^)9±(C;^)*-9±(C';c^) ^^^^^ 
= Jd'cJ d'C K{c,C;io)9±{Ccjr-g±{C';u;) 

= {g^{cj\K{u;)\g^{cj)) = -^T_j£^"l;g](^) > 0, 
which is the generalization of the paraxial result H91|l . 

If we take g{C', lo) — j (^; lo) to be the physical current density of interest, then the homogeneous approx- 
imation 

Xj{C,oj) = J d^C'D{CX';Lo)j^{C';io) (197) 

agrees exactly in its outgoing components with the actual radiation: 

XT'iC;Lo)^aiC,Lo), (198) 

but in order to actually satisfy the homogeneous Helmholtz equation everywhere, it must possess an equal 
magnitude of net power in the corresponding (i.e., conjugate) incoming modes making up X'"(C; w)i power 
which will flow in the opposite direction through any surface enclosing the sources. This extra power 
appears as extra work which would be done by the actual source charges on these incoming fields, if they 
were present. (These advanced fields do not lead to absorption by the sources making up j because, 
again, they are the time-reversal of the retarded fields which would be produced by —j* , not j itself.) 
Using Green's identities and the governing Helmholtz equations, one can show that this homogeneous 
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solution may also be formally expressed in the KirchofF-integral form, only with the Green function 
G'"*(C; C'; t^) replaced with the fundamental homogeneous solution D{C;C;uj) : 



XjiC;co)= J da'[a(c';c.)(n'.9')^(C;C';^)-^(C;C';^)(n'-a')a(C';^)]- 



(199) 



Unlike the inhomogeneous case, this remains valid for any boundary radius ^3 > and any observation 
position ^. 

Not surprisingly, the homogeneous approximation XjiC'j^) to a{(^;oj) is just the result of the above- 
mentioned one-to-one mapping between the isomorphic spaces Hout and Hhomi involving the replacement 
of each "outgoing" spherical Hankel function h^p [kQ with the corresponding "standing-wave" spherical 
Bessel function ji{kC) in the spherical wave expansion for a(^; lo), without altering any of the multipole 
expansion coefficients or basis vector fields for the angular dependence. Schematically, we have: 



Whom 

x7M - ci^lM. 



(200a) 
(200b) 

(200c) 

(200d) 
(200e) 
(200f) 



We expect that this X^{C'^ ^) is in some sense the closest possible homogeneous approximation to a{C,\ uo), 
which will be confirmed below by the non-paraxial form of the Maximum-Power Variational Principle. 

We can also express the various radiated power spectral densitities directly in terms of the multipole 
expansion coefficients or associated inner product, by using (|187|) in the limit C2 00. Asymptotically, 
i.e., for kr ^ £, the limiting forms of the spherical Bessel, Neumann, and Hankel functions are: 



neikr) 
hf\kr) 
hf\kr) 



sin(fcr- ^) 
kr ' 
cos(fcr - ^) 



kr 



^-\-ikr 

kr 

^ — ikr 

kr ' 



(201a) 
(201b) 
(201c) 
(201d) 



Examining the form of the Sommerfeld radiation conditions, we see that, asymptotically as k( ^ 00, 
any vector potential xp{^;uj) (not necessarily a source-free solution) may be locally decomposed into its 
incoming and outgoing components on the spherical surface i9-B(C) by the projections: 



^°"'(C;^) 
V^'"(C;^) 



k ec 



Equipped with these results, it is straightforward to show that for any solution tp E Haut 
outgoing power spectral density is given by 



(202a) 
(202b) 
H,„, the 



duj 



5'em[^°"1(oo;w) 



lUJ 



lim 

C2— >oo 2 



(fa [t/j""' • (n • a)i/)°"'* - 1/)°"'* • (n • d)ip'""] 



lim 

C2— >oo 



= lim 



9B{Q 
ILO 



J [(1 



k dC 



)V'-( 



dC k dC 



-I- c.c 



(203) 



dB{C2) 



^Cl/ J d9 sin4(l-t^)^.(^-t^)V'*] 
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which, after a httle algebra, reduces to the inner product over the outgoing multipole coefficients: 

,[^-](oo;c.) = [lCr'(^)|' + = ir^'i^) , > 0, (204) 



duj ' 



with equahty (vanishing outgoing power) in the last step if and only if i/j°''*((^; vanishes identically. 
Similarly, we find 



-^T™[^'1(oo;u;) = (^^''^c^) , V''°(c^)).„ > 0, (205) 



so we may write 



(xPicu) , xPiLu)) EE ir-'iLo) , V'°"'(c.))„, + (rP'" (u;) , V'-('^))„, 

= £ Tem [^"1 (oo; c^) - £?™ [•0-] (cx); cu) (206) 
= ||;5',m[^°"'](oo;c^)| + |£?™[^'1(oo;c^)| . 

The inner products previously introduced for Tiout a-nd 7ii„, and Tiom 7ii„ are just equal, respectively, 
to the (scaled) outgoing power spectral density, the absolute value of the incoming spectral density, and 
the magnitude of total (not net) in-flowing and out-flowing power spectral density, passing through an 
arbitrarily remote boundary somewhere outside the support of the (real or effective) sources. Choosing 
— Xji the relations H2()4|l and H2()5(l provide an explicit verfication of (|195|l . As an aside, we note that 
these separate incoming and outgoing power spectral densities and their (incoherent) sum, are often of 
more physical interest than the net power ^ J'em[''/'] (oo; w). For UV and lower frequencies, detectors 
(such as photomultipliers or solid state photoelectric devices) are typically directionally sensitive, in 
the sense of an acceptance solid angle J7 < 2tt, so measurements of the power in the radiation tp will 
typically reveal something closer to either |^5'EM['0°"'](oo;a;)| or |^Tem['0"'](oo;w)| , depending on 
the orientation, (or some fraction thereof, with less than full Att solid-angular coverage of the localized 
source), rather than ^^Pem [''/'] (oo; w). In the X-ray regime, shielding is difficult and detectors are often 
not directionally sensitive, and the response will likely depend on the total fluence into the detector 
volume, i.e., on \-£^ysM[4>°'"]{oo;uj)\ + |^J'em['/''°](oo; w)| rather than ^Tem['0](oo; w). 

In a similar fashion, we may formally express the spectral density of mechanical work in terms of the 
multipole coefficients for the exact radiation a and the homogeneous approximant x, by using H188|l . 
choosing a spherical boundary dB{(^2) and letting C,2 oo. After the algebraic dust settles, we find 

-l£y_,^[e^^;j]{u;) = Re^ [x,^,^(u;)*a,\,(c.) +xr™(^)*ar™(^)] 

em (207) 

= Re(x°"'(^), aM)„, 

In effect, we have managed to re-express the three-dimensional, "volume" inner product {j {uj)\a(Lu)) in 
terms of the two-dimensional, "boundary" inner product (x°"'('^) i '^(^))out ■ choose X = Xj, then 

by applying 1204(1 . the power spectrial density in the form ((207|l is seen to satisfy ((196(1 explicitly. 



E. Variational Principle 



Our general variational principle follows directly from the nonnegative-definiteness of the operator K{uj), 
i.e., equation l(173() . and the above energy conservation considerations. Re-tracing the standard proof of 
the Cauchy-Schwarz inequality, we find that, with a little extra care, it goes through almost unchanged if 
we replace the positive definite inner product with a positive-semidefinite pseudo-product. In particular, 
with the bilinear form associated with K{uj), one can show 

Mu;) \Kiu;)\g'{cj))f < {g{u;) \K{u)\g{u)) (g'(c.) \K{u)\g'{u)) , (208) 

with strict equality if and only if there exist constants (5, 5' , 5q e C, not all zero, and a non-radiating 
source go(C;'^) which is not identically zero but is otherwise arbitrary, such that 



5 g(C; a.) J' g'(C; c^) + <5o 9o(C; '^) = 0. 



(209) 
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We also have, for any complex number c € C, the inequalities 

Re[c] < |Re[c]| < |c| , (210) 

with overall equality if and only if c is real and nonnegative. Choosing g' = j± (^; ui) to be the solenoidal 

source responsible for the actual radiation a(^;cij) and g = g±{(^;Lo;a) to be an effective source cor- 
responding to a solenoidal, source-free trial vector field xiCl '^j oc) depending on certain variational 
parameters collected in the vector a, these inequalities yield: 



1/2 r a a. r„i//- .. .Ml/2 (211) 



with overall equality assuming ^ J'EM[a](C2; 7^ 0, if and only if 

X°"*(C;^^;a)cxa(C;a;). (212) 

(In the trivial case where -^'J'-em[0']{C2', 1^) = 0, equality is achieved trivially for xiC'i^) = or indeed 
for any x(C;'^) which is Hilbert space orthogonal to j±{C]^)-) 

In practice, we can obtain an actual variational approximation with one of two different but ultimately 
equivalent procedures, each arriving at the same final result when optimizing over equivalent families of 
trial solutions. In the first approach, we may partition the variational parameters as a = (%, ^O; ot'V ' 
where a' = (a'^, a'2,...,)^, and consider a restricted family of solenoidal, homogeneous trial solutions 
v{<^;u; a') normalized to unit power: 



:3'EMK'](C2;^;a') = 1, (213) 



or really to any conveniently-chosen but constant value. (Recall that this power is independent of the 
precise choice of surface F(C2) as long as (2 > Ci-) We then maximize |^ymech[eD_L; j](w; a')\ for fixed 
output power, in order to determine the optimal parameter values a = Q;'[j](a;) and therefore the 
relative shape (spatial profile and polarization) of the radiation. Once the relative profile is determined, 
we then renormalize this trial solution, taking 

X{C;^;a) = fioe'^'>viC;uj;a'), (214) 

for some real phase Oq, < 9q < 2tt, and nonnegative amplitude 770 > 0. To maximize the mechanical 
work term, the overlap integral (3D inner product) between (scaled) solenoidal electric field and current 
density must necessarily be real and negative, so ^0 is chosen so as to ensure 



■iRe 



Im 



{e^^{Lu;cx')\jiL,))^>0 (215a) 
{s^^{u);dc'\j{u))]=0. (215b) 



Finally , the positive amplitude 770 > is chosen to ensure energy balance between the power (spectral 
density) in the approximate outgoing fields and the work which would be done on these fields, if present 
in the region of the sources, by those actual sources: 



-ir,oe-^^° {s^^{u,)\j{w)) = 4£?^u[v°''^]{(2;u,;a), (216) 



or 



~ _ 2 C)lJ- mccnL-^-i-,^jv-, — / (217) 

£5'BMM(C2;a;;d) 



Note that under the assumptions of energy conservation given the sources, the overall amplitude, set by 
f]o, is uniquely determined, and the overall phase is determined modulo 27r, because the mechanical 
work done on/by the sources is linear in the fields, but the radiated power is quadratic in the fields. 
As in the paraxial case, we can also interpret this as a maximization of the small-signal, no-recoil gain 
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which would be experienced by the trial field if it were to be incident on the sources. Upon substituting 
back into the Cauchy-Schwarz inequality, canceling common terms, and squaring, we deduce 

-i£5'_„[£^^;j](^;a) < ^T™[x°"'](C2; c^; a), (218) 
with equality if and only if 

x™'(C;w;a) = a(C;^), (219) 

or equivalently, 

x{C,cj;cx)^Xj{Coj)= I d^C D{C;C-:^)j±{C:U^)- (220) 

Thus, we sec that the variational approximation x°"'(C2; oj; a) may also be obtained more directly from 
a simultaneous optimization, with respect to all parameters in a, of the outgoing power, subject to the 
constraint of energy conservation: 



a = argmax 



|TyEM[x°^(C2;^;a) (221a) 



s.t. -i£^a'_j£^^;j](a;;a) = £^J',„[x°"'](C2;^;a). (221b) 

This is our desired result. As anticipated, we have arrived at a variational principle which holds in the 
general, non-paraxial case, which is exactly analogous to that found previously in the paraxial limit. 
Unfortunately, there seems to be no simple way to express the constrained optimization problem (|221|l 
as an unconstrained maximization of a ratio between the mechanical work and radiated power (or some 
functions thereof). 

It is also instructive to also see how this all works out explicitly in the spherical wave basis. Applying 
the usual Cauchy-Schwarz inequality directly in Haut, and using equations 1)204(1 and 1207(1 . we have 

-i£j'_j£^i;j](c.;a) =Re(x°"'(^;a), a(^))„„^ < |(x°"'(^; «) , aMLJ 

< ix-^co; a) , x°"'(^; («(^) , (222) 

= IjJ'bm (^; c^)'/' £j J'bm [a] (ex.; l,)'/^, 

with overall equality if and only if x°"*('^;ck) = k{ljj) a{uj) for some positive real scalar k{uj) > 0, or 
equivalently when 



> 0. (223) 



for all mutlipole coefficients of a{(^;uj). Imposing the additional energy conservation relation on 
X((J;w;a), this becomes 

-i£j'.„ecj£xi;j](^;«) < 5'EM[a](oo;c^) (224) 

with equality if and only if x°"'(i[';w; a) = a(^; w). Because in practice we will not have knowledge of 
the actual multipole expansion coefficients for a in situations where we are resorting to a variational 
approximation, this "angular-representation" or "boundary-value" formulation is not directly applicable 
for calculation, but it does lead directly to the relevant power inequality without any worries about a 
semidefinite or "pseudo" inner product. We also see that the MP VP follows immediately as a consequence 
of the minimum distance property of orthogonal projectors in these multipolar Hilbert spaces: the 
variational approximation xiCl is simply that which minimizes the distance to the actual solution: 



a = argmm 



(a(c^) - x{^; a) , a{Lj) - x{^; cx)f^] (225) 
within the considered family of solenoidal, source-free, trial wavefields. 
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V. DISCUSSION 

Here we assess advantages and limitations of the MPVP, briefly consider possible variations and exten- 
sions, and compare the method to other variational techniques commonly encountered in electromagnetic 
theory, in particular those which have been used for paraxial radiation problems, laser-plasma interac- 
tions, and FEL-type devices. 



A. Summary 

The derived MPVP H221|l . valid in the paraxial or non- paraxial regimes, says, in effect, that classical 
(accelerating) charges spontaneously radiate "as much as possible," consistent with energy conserva- 
tion. By spontaneous we mean that charges are assumed to follow prescribed classical trajectories 
determined by the external fields and possibly (averaged) quasi-static self-fields, but are not affected by 
recoil/radiation-reaction, multiple scattering, absorption, or other feedback from the emitted radiation 
(although arbitrary pre-bunching due to previously-imposed or radiated fields may be included, if it can 
be characterized). The sources are assumed to be localized in space, so that the far-field (i.e., radiation 
zone , or wave zone) is defined, and at least weakly localized in time, so that Fourier transforms (i.e., 
frequency-domain representations) of the current density exist. In turn, the radiation, once emitted 
by any part of the source, propagates as in free space, without further scattering or absorption. The 
optimized trial mode-shape (or more accurately, the outgoing component thereof) is then the best guess, 
within the manifold of possibilities allowed by the adjustable parameters (a), of the exact radiation field 
profile radiated by those prescribed sources, and the optimized outgoing Poynting fiux yields a lower 
bound on the actual spontaneously-radiated power spectral density at the frequency under consideration 
(apart from any numerical/roundoff errors in the computation). The approximation will improve mono- 
tonically as additional independent parameters are included to allow for more general envelope shapes. 
As with other variational principles, the relative error in the value of the variational functional (here, 
the power) at the stationary point is generally smaller than the error in the parameterized trial function 
(here, the radiation profile). The adjustable parameters may appear linearly or non-linearly in the trial- 
functions; if they are taken to be the linear expansion coefficients in a sum over orthogonal modes, then 
the MPVP is equivalent to a truncated basis-set approach. (In the paraxial case, "outgoing" implies 
propagation nearly along the optic axis away from the localized sources, defining the "downstream" or 
"post-source" direction, while "incoming" means propagating towards the sources from the upstream 
direction. In this case, say, a right-going pulse has both ingoing and outgoing components. In the gen- 
eral non-paraxial case, "outgoing" implies wavefronts which, at least asymptotically, diverge in time and 
propagate away from the localized sources, while "incoming" denotes the reverse, i.e., asymptotically 
converging wavefronts propagating towards the sources.) 

This variational principle can be variously interpreted according to one's tastes or application. As 
we have seen, the best variational approximation maximizes the radiated power consistent with the 
constraint that this energy could have arisen from work extracted from the actual sources. It also 
minimizes a Hilbert-space distance between the actual fields and the parameterized family of solenoidal, 
free-space radiation fields, and in many cases may actually be regarded as an orthogonal projection into 
this manifold of trial solutions. It also maximizes, for each frequency component, the spatial overlap, 
or physical resemblance, between the actual current density and the trial fields, when the latter are 
extrapolated back into the region of the sources assuming source-free propagation. Equivalently, one can 
say the optimal free-space field profile is that which, if it actually were to be incident on the sources, 
would maximally couple to the given sources and would experience maximal small-signal gain, and 
furthermore the virtual gain delivered would be equal to the estimated power spontaneously radiated. 

Whatever the shape of the trial function, and whether the variational parameters appear linearly or 
non-linearly, after optimization this profile may be regarded as one of a complete set of solenoidal basis 
functions solving the source-free wave equation. After introducing an inner product associated with 
the radiated power (not energy), these basis functions can be orthonormalized via a Gram-Schmidt 
procedure, leaving the original trial function unchanged. Bessel's inequality then confirms that the 
outgoing power radiated in (the outgoing component of) this one mode can be no greater than the total 
outgoing power radiated in all the modes. All this is obvious. What is perhaps more surprising (if we do 
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not think too hard about the Fundamental Theorem of Calculus or its multidimensional generalizations) 
is that we can exactly relate the outgoing power radiated in the far-field to the "virtual" work which 
would be performed by the localized sources on the "source-free extrapolant" of this radiation field, as 
if it and it alone were present in the region of the source, i.e., to the overlap integral between source 
and extrapolated electric field, even though the extrapolated field may not closely resemble the actual 
emitted field there, as it must necessarily solve the inhomogeneous wave equation. The factor of i 
in the resulting conservation constraint emerges in order to avoid double-counting in the energetics. 
The radiated power in the outgoing component of the variational approximation is being related to the 
mechanical power which would be delivered by the sources to the full source-free fields, if they were 
actually present in the region of the sources, rather than the actual, inhomogeneous fields; and, any 
wave-field satisfying the source-free wave equation everywhere in space must have the power in each 
outgoing mode exactly balanced by that in a corresponding incoming mode, or else singularities would 
arise somewhere in space. (Without any sources, in the paraxial case, radiation observed to emerge 
from, say, the right side of any longitudinal interval must have previously traveled into that interval 
from the left. In the free-space non-paraxial case, analogous conclusions can be drawn for plane waves 
traveling along any chosen direction, or, thinking in terms of multipole radiation, it is necessary that 
any diverging spherical wave must be matched to a converging spherical wave to avoid a singularity at 
their common center of wavefront curvature.) 



B. Assessment 

Again, the present approach involves approximation of the actual radiation fields by trial modes which 
are solenoidal and satisfy the homogeneous (source-free) wave equation. As we are primarily interested 
in the net radiation from the given sources observed in the free-space region beyond the sources, it is 
appropriate to use solutions of the source-free equation there, but then it is natural, efficient, and even 
almost inevitable, in the context of an approximation scheme, to effectively extrapolate the fields back 
into the support of the sources by free-space propagation. (Otherwise, if we could somehow propagate 
according to the actual dynamics, we would not need to resort to any approximation - but we would also 
be carrying around unneeded information about the near- fields.) In the paraxial limit, these free-space 
solutions are uniquely specified by the carrier frequency and the (complex) profile in any transverse 
plane, which can be decomposed into a convenient, countable set of modes consisting of the eigenmodes 
of some quantum-like operator(s). In the general case the solutions are more complicated, but from the 
spherical wave expansion, we know at least that they also form a separable Hilbert space, parameterized 
by the multipole expansion coefficients. Because there is no known closed-form, analytic solution for 
a general focused radiation beam outside the paraxial approximation, the formalism and variational 
principle are probably of most use in the paraxial case, although may also be applied, with little added 
complication, to cases where the radiation consists of a superposition of multiple beams, each paraxial 
individually, but with sufficiently divergent directions to violate overall paraxiality, or to cases where 
higher-order terms in the generalized paraxial (i.e., diffraction-angle) expansion are included to describe 
a moderately-collimated beam. 

The output x(^;a;; a.) of the MPVP is therefore actually better considered a variational approximation 
to the source-free "embedding" or "extrapolation" XjiC'j'^) rather than to the actual outgoing vector 
potential a{C;uj). If the radiation is sufficiently directional, then the outgoing component x°"'(^; a) 
may be determined simply by restricting attention to some finite solid angle. In the most general 
case, where the "spurious" incoming radiation may overlap with the "physical" outgoing radiation, it is 
more difficult to disentangle and x'" ^-t an arbitrary position. If the spherical wave expansion is 
known, then the outgoing projection may be easily determined everywhere in space, but direct use of 
the the spherical- waves as the variational basis is not generically expected to prove convenient. (This 
would correspond to just using a truncated multipole expansion for the radiation.) In the far-field 
limit fc^ — > oo, at least, the outgoing component may be determined easily by using the Sommerfeld 
projections (|202|l . 

Because of these difficulties, and remaining mindful of the proportionality between the calculated power 
delivered from the source J_l(C;^) to the radiation fields corresponding to the vector potential a{(^;uj) 
or to the fields of its source-free approximant x(^,lj), one might be tempted to apply the variational 
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method directly (just without the extra factor of i) to outgoing, solenoidal solutions of some conve- 
nient inhomogeneous or homogeneous Helmholtz equation. However, purely outgoing-wave, free-space 
solutions cannot exist everjrwhere in space, and will inevitably possess one or more singularities in the 
vicinity of the sources, preventing calculation of the required mechanical work integrals. Furthermore, 
any solenoidal vector field whatsoever, say is trivially the solution to some inhomogeneous 

Helmholtz equation for some boundary conditions and some solenoidal source, namely, for those bound- 
ary conditions satisfied by itself, and for a solenoidal source given by —(9^ + )'(/'. Without further 
constraints, this class of trial solutions is simply too general to be useful. If we naively apply the MPVP 
using homogeneous trial functions, the method will converge not to a but to something proportional 
to j, since, for given norm, no other vector field can have greater overlap with j than a multiple of j 
itself. By using the source-free solutions, we effectively regularize the problem, so that we can converge 
to something which maximizes the overlap with the source under the assumption that the trial look like 
a radiation field, and not a current density. 

In any case, the method requires not only the parameterization of trial modes but the calculation of 
inner products between these trial modes and the current density, or else some other, problem-specific, 
means to determine the absolute level of the power, either analytically, in terms of the free parameters, 
or numerically and indeed repeatedly, as the parameters are varied in search of the maximum. One 
can easily imagine counter-examples where it is simply easier to work directly with the Lienard-Wichart 
potentials or corresponding expressions for the fields, than use the variational principle. 

We also have some options (or rather, trade-offs) in handling stochastic sources. So far, we have assumed 
a specific, prescribed, deterministic current density t), but in practice we must contend with statis- 
tical uncertainty in the exact electron trajectories in any one realization, meaning ji^Ci'-^) is in principle 
stochastic, with statistical properties determined by the phase space distribution of the particle beam (in 
the reduced six-dimensional phase space, if we may neglecting direct two-body or higher-order interac- 
tions, but in an even higher-dimensional phase space if particle-particle correlations are to be included). 
In certain cases, we can get by with using an averaged or course- grained current directly as the source, 
without need for further probabilistic considerations. In the Vlasov equation, for example, the fields that 
appear arc self-consistent mean fields, i.e., solutions to Maxwell's equations for the averaged sources. 
The exact fields arc linear in the sources, so it makes no difference in the end result whether we average 
the currents first and then calculate the resulting fields, or calculate the fields for a general instance of 
the currents and then average the result - the averaging procedure will commute with the convolution 
over the Green function. (Pragmatically, the former approach will almost always be more efficient if the 
average is all we need, since the random variables determining the stochastic shape of the current need 
not carried through the convolution, but the latter approach allows calculation of higher-order correla- 
tions or moments.) When using the MPVP to approximate the radiation, the averaged trial fields will 
equal the trial fields from the averaged sources if the variational parameters appear as linear expansion 
coefficients in a basis-set expansion, but not, in general, if they occur non-linearly in the trial wavefields. 
Convolution is a linear process, but optimization is in general a nonlinear one. 

Unfortunately, it is often the case that we desire to estimate the fields from given sources, then calculate 
moments, rather than to determine the fields from averaged sources. This is particularly true of spon- 
taneous wiggler radiation, where the radiation is almost entirely "fiuctuational," both in the sense that 
the squared-magnitude of the average field is small compared to the average of the squared-magnitude, 
and in that the relative variance in the spectral intensity approaches 100% in each sufficiently narrow 
frequency-band for radiation from sufficiently long, unbunched electron beams. At least second-order 
moments are needed to estimate the power spectrum and transverse coherence properties, and fourth- 
order moments are necessary if we wish to to estimate the uncertainty in the estimates for these spectra 
or coherence functions. Statistical estimates for the power spectral density, based on the moments of 
the MPVP trial fields, cannot in general be guaranteed to be lower bounds, as in the deterministic case. 

Before even calculating moments, we may cither approximate the radiation fields from the total current 
density using the MPVP, or else apply the MPVP separately to each particle trajectory or family of 
trajectories, then superimpose the resulting variational fields. If identical, linear basis-sets are used 
in the expansions for each particle's radiation, identical results will be achieved (although the amount 
of calculation or computational effort may differ), but if nonlinear variational parameters are used or 
even if different truncated linear basis sets are used for different particles, the two approaches will be 
inequivalent. 
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A simple example should clarify the issues. Suppose we seek to approximate the wiggler radiation 
from some electron beam with negligible normalized energy spread ^ 1 but some non-zero transverse 

emittance e > Ao, where Aq ~ 2j^'^'>^^^~^ I'^u) central wavelength of the wiggler radiation in the lab 

frame, A^ is the wiggler wavelength, and a„ is the normalized wiggler strength parameter. For radiation 
from a single on- axis electron, confined within the relative bandwidth ^ ~ where -/V„ is the number 
of wiggler periods, we know the radiation is approximately diffraction-limited, with characteristic angle 
SOq ~ ^ and a waist located at the midpoint of the wiggler. However, the radiation for the beam as 

a whole will not be diffraction-limited, but will have a transverse degree of coherence r<,„i, j. ~ ^ < 1, 
due to averaging, or convolution, over the transverse particle distribution. 

Suppose first we use the total current density as source, and use paraxial Gauss-Hermite or Gauss- 
Laguerre modes as trial solutions, with optic axis coinciding with the wiggler axis. If we attempt 
to include, for A « Aq, only the fundamental Gaussian mode with undetermined waist size wq and 
location then, regardless of the exact values of the optimized parameters, the approximation cannot 
simultaneously predict both the transverse size and angular divergence with accuracy, or equivalently, 
the partial transverse incoherence cannot be adequately captured, since we are putting all the power in 
a single, diffraction-limited transverse mode. In fact, we expect that we must include at least n_L ~ 
diffraction-limited modes to resolve the transverse coherence, even if the average transverse intensity 
profile continues to look Gaussian. If the particle phase space distribution function is also taken to be 
Gaussian, then the needed overlap integrals in the MPVP may be calculated analytically. 

Alternatively, we could decompose the current density into contributions from each electron, and employ 
a single Gaussian mode for each, only with the optic axis aligned along the average single-particle 
trajectory determined by its spatial and angular displacement from the reference orbit. The needed 
overlap integrals in this case are easily perfomed, at least for typical values of Aq, 7, and a,j, where 
the transverse spatial extent of an electron orbit is small compared to wq, but after determining the 
single-particle MPVP fields, each with a different propagation axis, they all must be superimposed to 
determine the full field. This can be done easily in the far-field so as to determine the angular spectrum, 
or at arbitrary positions, to a good approximation, provided the angular deviations for the electrons 
remain moderately small, i.e., < i. The result will agree approximately, but not exactly, with the 
previous approach. The details of a simplified version of this calculation will be presented in elsewhere. 



C. Extensions 

The results derived here are restricted to radiation from localized sources otherwise propagating in 
vacuum, which includes the cases of charged particle beams in bending magnets or undulators, but it 
would be desirable to extend applicability to more general radiation problems. Using the macroscopic 
Maxwell's equations, it should be straightforward to include a linear, dielectric background eoCw) > 
which is spatially homogeneous, but possibly frequency-dependent. This would be sufficient to handle 
many Cerenkov radiation problems. To treat some forms of transition radiation and waveguide prob- 
lems, piecewise-slowly-varying spatial variations in the dielectric tensor and the possibility of (perfectly) 
conducting boundaries should both be included, but in general the "source-free" solutions become dif- 
ficult to calculate and parameterize - a single planar jump, or other simple geometries, might still be 
tractable. 

One important class of geometries for which this formalism should be almost immediately applicable is 
that of conducting-wall and/or dielectric waveguides with translational symmetry along the optic axis. 
In such circumstances, it is well known [ll. Isl lsS. l40l | that such systems share much of the mathematical 
Hilbert-space structure as that elaborated above for paraxial radiation - in particular, solenoidal modes 
which are uniquely specified by their cross-section in any transverse plane, and which are orthonormal- 
izable with respect to the natural £<2/Euclidean inner product in that plane. (However, these modes 
are not complete for the near-fields, as discussed, for example, in addition, it is easily shown 

that, in any such waveguides, the power in any source-free, propagating mode, measured sufficiently far 
downstream from a localized source inserted inside the waveguide, can be related to the overlap integral 
between the actual current density of the source and the electric field profile of this free-space mode 
extrapolated back into the support region of the sources as if no sources were present. Modifications 
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due to finite conductivity efi'ects might be then treated perturbatively. Again, this is all familiar, and is 
exactly analogous to our basis-set expansion technique. What is perhaps not widely appreciated is that 
the Cauchy-Schwarz inequality then allows one to transform this basis-set approach into a maximum 
principle, where the exact form of the normal modes are not needed. 

Attempts to use a more general complex- valued dielectric tensor or index of refraction in order to incor- 
porate linear loss or gain will encounter fundamental difficulties. Our derivations have exploited the fact 
that the free-space retarded and advanced Green functions are just related by complex conjugation in the 
frequency domain, but with any dissipation or gain, this time-reversal symmetry will be lost. Ideally, one 
would like to fully generalize the treatment to allow for dielectric tensors e^^(a;, t; lo) with any frequency 
dispersion consistent with the Kramers-Kroning relations and/or piecewise-slowly- varying spacetime de- 
pendence (inhomogeneity) , and possibly anisotropy, in order to model the effects of essentially arbitrary 
conducting or dielectric boundaries, obstacles, or passive optical elements, such as lenses, prisms, polar- 
izers, waveplates, apertures, waveguides, cavities, transmission lines, fibers, gratings, etc., either ideal or 
lossy, but it not even clear whether this is possible, or if so, practical. Many variational techniques have 
been used for waveguides, cavities, and other structures, but these are usually of a stationary, rather 
than extremal, character, an important distinction which is examined below in Section [V El All of these 
directions are left as open problems. 



One is also naturally led to speculations as to whether similar ideas may be applied to stimulated 
emission. Naively, one might imagine some iterative procedure, where the approximated radiation, 
assuming given sources, is somehow used to deduce corresponding energy changes that must have taken 
place in those sources, which leads to revised estimates for the radiation, and so on. It is not clear 
whether this is applicable to any real problems, or even if it can converge to a self-consistent solution 
describing an active medium. 

However, results for the case of spontaneous emission from a classical particle beam may be of direct 
relevance to the problem of stimulated emission, namely in the so-called small-signal regime where 
saturation and depletion effects are ignored. In his celebrated 1917 quantum analysis of radiation j4^ . 
Einstein first classified the radiative processes involving photons (or really any bosons) into spontaneous 
emission of, stimulated emission by, and stimulated absorption of the photons by atoms or other material 
sources. By exploiting the fact that the matter and radiation can be in thermodynamic equilibrium, 
Einstein used the properties of blackbody radiation and the principle of detailed balance to establish the 
proportionalities between the intrinsic rates, within each mode, for each of these processes, as summarized 
in the famous A and B coefficients, results which must hold for arbitrary initial conditions, not just those 
consistent with thermal equilibrium. Less widely appreciated is that these relations persist in a purely 
classical limit of the matter and fields where Planck's constant h never appears, clearly discussed, for 
example, by Beckefi in Chapter 2 of ji^l, and derived using somewhat different arguments in [il ]. 

As a simple example, consider the case of radiation transport in an electron plasma, where direct many- 
body effects are neglected, and with an isotropic, single-particle momentum distribution function f{p), 
where p = \\p\\ is the magnitude of kinetic momentum, and uniform density uq over some region. (This 
might also describe a long electron beam in its average rest frame.) Then, for each type of radiative 
process (or the aggregate of all processes) , the spontaneous emission rate i^i defined as the power emitted 
in polarization e per unit volume of medium per unit bandwidth per solid angle in the direction of s, 
may be written as 



where r]^{p] e; s) is the single-particle, intrinsic emission coefficient, which is determined by the details of 
the microscopic physics, but is independent of the particle distribution function or the incident radiation 
state. The net absorption rate, i.e., bare stimulated absorption less stimulated emission, is given by 
aojluj, where the radiant brightness is the incident power per unit area per unit solid angle per unit 



D. Connections to Stimulated Emission 




(226) 
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bandwidth, and the absorption coefficient is given by 

a^[e;s]^ — T-^no / d^pT]^{p;e;s)-§^f{p), (227) 



where e — e{p) = y^?p^ + m^c? is the particle energy, rij, is the background index of refraction, and 
rj^oip] e; J) is the same quantity that appears in the spontaneous emission. These results are completely 
classical, as no factors of h anywhere appear. The dependence on the derivative of the distribution 
function is the only reminder that this expression represents the net difference between stimulated 
absorption and emission. These connections between spontaneous and stimulated emission/absorption 
are essentially a generalized manifestation of the well-known fluctuation-dissipation theorem, which 
relates the response of a system when perturbed to the spontaneous fluctuations which occur in the 
absence of external perturbation jiM \aA l45l | . 

In FEL physics, first described quantum- mechanically but now understood classically, such generaliza- 
tions of Einstein's arguments are the basis for Madey's theorem in one-dimensional theory |46| and its 
generalizations to higher dimensions (see, for example, ^3,EHE3i or else ^3 references therein), 
where the gain curve (specifically, the relative change in intensity versus de-tuning) in the small-signal 
regime is proportional to the derivative of the spontaneous emission spectrum. 

Given these connections, certain properties of the spontaneous wiggler radiation, or approximations 
thereof, can yield information about the stimulated emission, or gain, in the small-signal regime. Perhaps 
more satisfying is the obverse relation, whereby our intuitions about stimulated emission provide the 
clearest justification and interpretation for the appearance of a maximum-power variational principle for 
classical spontaneous emission. In situations with gain, we might naturally approximate the radiation 
profile by that which maximizes the extraction of energy from the electrons to the laser. For example, 
this is precisely the heuristic justification for the variational principle suggested in for optically- 
guided modes in FELs, where the fundamental mode is estimated by maximizing the imaginary part of 
the effective wavenumber. If the radiation is to be estimated by one mode, the mode should be chosen 
to be whatever shape is expected to have the fastest growth (or slowest loss). 

In fact, we saw above that the MP VP can be interpreted precisely in this manner - as finding the 
mode shape which, if actually incident, would maximally couple to the given sources, and furthermore 
the "virtual" gain delivered would be equal to the estimated power spontaneously radiated. The only 
essential difference between our case, and say, Madey's theorem is that by taking completely prescribed 
sources, we implicitly assume that any radiation, once emitted by some part of the source, can then 
propagate to the far-field where it may interfere constructively or destructively with other radiation, but 
cannot cause any recoil in that source nor subsequently be re-scattered or absorbed by any other part 
of the source. So in fact we find a relationship between the spontaneous emission spectrum and that of 
the hare stimulated emission, not the net response given by the difference between stimulated emission 
and absorption as in Madey's theorem. 



E. Relation to Lagrangian Formulation and Other Variational Approaches 

A survey of variational techniques in electrostatics, electrodynamics, and optics reveals that, at some 
fundamental level, they seem to derive from only a very few general concepts: Hamilton's principle for 
systems derivable from a Lagrangian/action, energy conservation/optimization, electromagnetic reci- 
procity, and entropy maximization (for problems of a thermodynamic nature.) Of course these principles 
are interrelated ~ both energy conservation and reciprocity follow from the structure of the governing 
Lagrangian, while entropy maximization usually relies on energy conservation or other dynamical in- 
variants as constraints and automatically satisfies certain reciprocity relations - but nevertheless they 
remain useful and conceptually distinct organizational categories. Each variational principle may be 
further classified by whether it demands (generally constrained) optimality, or mere stationarity, of the 
relevant function or functional. 

We examine in turn the relevance or relation of each of these principles to the MPVP. Obviously the 
MP VP is most closely related to energy principles, although in the case of classical spontaneous radiation, 
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we maximize the rate of energy transfer (power) in each spectral band, rather than, say, minimize 
potential energy as in electrostatics, or free energy in thermodynamics. 

Entropy maximization might be a useful technique to determine or approximate certain source distri- 
butions given only macroscopic thermodynamic constraints or other incomplete information, but is not 
directly relevant for the problem of finding the radiation from given deterministic sources, which is more 
a problem of dynamics than thermodynamics. However, as we saw above, new connections and rationales 
are revealed by generalization of the arguments leading to Einstein's A and B coefficients to classical 
beams or plasmas, and these arguments rely on the equilibrium properties of black-body radiation, even 
when the actual system of interest may be far from thermodynamic equillibrium. 

Reciprocity-based variational principles are commonly used to approximate resonant frequencies in cav- 
ities or waveguides, impedances of cavities, transmission lines, apertures, or other structures, and scat- 
tering cross sections by conductors or dielectrics 0, |^ 113 ■ They can be suitably generalized to many 
situations where the background medium is not necessarily homogeneous or isotropic. Despite superifi- 
cial similarities between the MPVP and reciprocity-based techniques, a closer look reveals them to be 
quite distinct. Consider the solenoidal electric field e^±{C;uj) derived from some Coulomb-gauge vector 
potential 'ip{C]u;), and any source j'(^;a;), not necessarily that associated with xp. For our purposes, we 
may define the complex reaction of the field £a_L on the source j' as 

n[e^^,j'] = J rf'\£^i(C;cc>).j'(C;c^) = (ev^J-MI/M) , (228) 

Which is similar to the original definition in 50], but involves only the solenoidal electric fields, and drops 
the analogous terms involving the magnetic field and magnetization density, because we are working 
with the microscopic fields from classical sources. Given any source j{C;i^), let a{C;u!) denote the 
corresponding Coulomb-gauge, retarded vector potential, and let xiC'-,^) ~ XjiC'-,^) denote the closest 
homogeneous approximant (i.e., that vector potential obtained by using j± as an effective source in 
convolution with the fundamental solution, rather than an actual source in convolution with the retarded 
Green function); and define a' and x' in analogous fashion for the source ji'. 

By using the symmetry properties of the Green functions, it is then straightforward to establish the 
Raleigh-Carson reciprocity relation for the inhomogeneous case: 

n[ea±,j']H = {ea±{ujr\j'{uj)) = {ea'±{ujr\j{uj))=n[ea'±,j]{io), (229) 

and similarly for the homogeneous case: 

n[e^^,j']{u;) = (e^^(a>)*|/(a>)) = {e^, ^{u;r\j{u;)} ^ TZ[e^, ^,j]{cj). (230) 

Such a quantity was introduced in (50| as a measure of "reaction," or coupling, between the sources j and 
j' , and in many contexts a variety of physical observables such as fields or forces, impedances, etc., may 
be proportional to or otherwise related to reactions. Despite resemblances, the complex self-reaction 
{e.^ ±{uj)*\j{u!)) is distinct from the complex power {e-if, ±{uj)\j{uj)) , because no complex conjugates 
appear in the overlap integral in the former case. (Adding to the potential for confusion, the imaginary 
part of the complex power is referred to as the "reactive power," following the conventions of circuit 
theory.) Said another way, the complex power is associated with a standard complex inner product, 
which is conjugate-symmetric and positive-definite, while the reaction is associated with a symmetric, 
bilinear form which is not a true inner product. Physically, the concepts are also distinct: energy 
conservation arises from the invariance properties of the full electromagnetic Lagrangian (including self- 
consistent particle dynamics) under time-translations, while reciprocity arises from invariance under 
time reversals. The real part of the complex power is of course proportional to the time-average (over 
an optical period) of the power exchanged between source and field. The real part of the reaction is 
actually proportional to the fluctuation about this average. 

The concept of reaction is of interest precisely because of reciprocity properties like H229|l or H230|l . 
and the possibilities for variational approximation that emerge, by demanding that all trial sources 
"look" the same to themselves as to the correct sources. For example, suppose we seek a trial-function 
approximation xp to the unknown vector potential xp, given the actual source j. Let j correspond to a 
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source associated with a (the existence of non-radiating sources impHcs non-uniqueness, which is of no 
consequence here.) From a reaction point-of-view, we would hke an approximation such that 

n[s^^,j]{w)^n[e^^,3]{u), (231) 

but by assumption wc have no tractable way of calculating the right-hand side, so instead that approx- 
imation is then "best" if we can impose the constraint 

^[£^±>j1H«^[£^±,i]H- (232) 

In fact, one can show that under this constraint, the reaction TZ\£^ ^,j\{u)) is then stationary for first- 
order variations of about the actual solution xj^. 

Note that this leads to a stationarity condition, not an optimization condition. Because reaction and 
power are not equivalent, in general they lead to different variational principles. In the Method of 
Moments, which encompasses both the Finite-Element and the Ritz-Galerkin basis-expansion tec;hniques, 
one can obtain algebraic (typically linear) stationarity conditions using either a symmetric form between 
trial functions and weight functions, and thus automatically satisfy reciprocity relations, or using a 
conjugate- symmetric inner product, and satisfy energy constraints, but typically not both simultaneously. 
Perhaps the reaction approach may also lead to some useful variational principle for our problem, but 
we do not pursue this question further here. One immediate difficulty is that, if we use inhomogeneous 
solutions ip as trial vector fields, the corresponding sources are easily determined by substitution into the 
Helmholtz equation, but the fields themselves are difficult to parameterize in any economical manner, 
while if we use source- free solutions for tp, variational families may be more easily characterized, but the 
effective sources appearing in the reactions are then difficult to determine. 

Now we turn to action principles. Because the fundamental dynamics of charge-carrying matter and 
electromagnetic fields are ultimately derivable from a Lagrangian, we had anticipated that our maximum- 
power variational principle would be traced back to Hamilton's Principle, but this appears not to be 
the case. (Perhaps this should not have been so surprising; after all, the Lagrangian involves differences 
between kinetic and potential energies, while an an energy-based principle should involve the sum.) To 
understand the issues encountered, let us proceed by rc-framing the governing dynamical equations in 
terms of a standard Lagrangian formulation. For prescribed current and charge densities, the dynamics 
of the fields are derivable from the action functional, which in terms of the Coulomb-gauge potentials 
and scaled coordinates may be written as a space-time integral over the Lagrangian densities for the free 
electromagnetic fields and the source-field interaction: 

S = JdrJcl^C [A+A„t], (233) 

where 

Co = ||£(C;T)f - ||6(C;r)f = ||-|:a(C;T) - d^rW " 11^ x a(C;r)f , (234) 

and 



= j(C; r) • a(C; r) - m(C; rmC, r), (235) 

subject as always to the solenoidal gauge constraint 

a-a(C;T)= 0. (236) 

Using the reality of the physical potentials and sources and the (Hilbert-space) orthogonality between 
the transverse vector field a(^; r) and longitudinal electric field — 9^(C; r) when integrated over all space, 
we can re- write this as 

S = jdrJd'C [(|:«*-|ra) + (a<^*-a<^) - idxa*)-{dxa) + {a*-j) - (<^»] , (237) 

where we may use j or j± interchangeably within the action functional because the current density only 

appears in an inner product with the solenoidal vector potential a. The equations of motion (Euler- 
Lagrange equations) for the vector and scalar potentials (and their complex conjugates) are obtained as 
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the conditions for this action remaining stationary under independent, infinitesimal variations (5a(^;T), 
da{(^; r)*, (5(/)(C; t), and S(f){C; t)*, which are fixed at the space-time boundaries (here taken to be at infin- 
ity), are consistent with the gauge constraint, and which do not destroy integrability of the Lagrangian 
density, but are otherwise arbitrary. 

By interchanging the order of temporal and spatial integrations and using the Parseval-Plancherel Iden- 
tity, we may express this action in the (scaled) frequency domain: 

S = Jd^Jd'^C [{iuay-iiioa) + {d(t)*-d(l)) - {dxa*)-{dxa) + {a*-j) - (?!>»] , (238) 

with additional gauge constraint 

a-a(C;w) = 0, (239) 

where we have used integration by parts on the time-derivatives, and where each potential and source 
term appearing in the action integral is now interpreted as the Fourier transform in scaled time (assumed 
to exist) of the function of the same name appearing in the time-domain version above, and where 
coordinate dependence has been suppressed for the sake of brevity and readability in the equations. 
Using standard vector identities, this may be written in the more transparent form: 



S = jdujd^da*- (^2 + w^)a + a*-j - cj)* (9^^) - (j)* + d-{(jfd(l) - a* xdxa)] 



(240) 



Demanding stationarity under independent, infinitesimal variations 5a{C,]uo), 5a{C,\u)* , (5(/)(C;i^), and 
5(j){C,'uj) satisfying the gauge constraint and fixed boundary conditions (as ||^|| oo and \uj\ oo) leads 
to the frequency-domain version of Poisson's equation and the wave equation in the Coulomb gauge. Now, 
by Gauss's theorem, the final term involving the pure divergence can be expressed as a surface integral 
over the spatial boundary (at infinity) , which will not affect the stationarity conditions (Eulcr-Lagrange 
equations) under the assumption of infinitesimal variations with fixed boundary conditions as dictated 
by Hamilton's principle, so it may be dropped without any change in the resulting dynamical equations. 
Likewise, we may then introduce the divergence d ■ q oi some other, arbitrary differentiable vector field 
q which may be a function of C, and lj and a functional of a and a* and their spatial derivatives. The 
significance of this auxiliary field will emerge shortly. Since we are ultimately only interested in the 
radiative component of the fields, the terms in the action involving the scalar potential (p may also be 
suppressed, as they are constant with respect to variations in the vector potential. Recalling that the 
time-domain potentials and sources are all real, we can then use as the relevant action the modified 
functional 



oo 

S = jd^jd^C [a*-{d'^ 



-uj'^)a + a*-j + d-q\+c.c., (241) 



and, noting the additivity with respect to contributions from the various frequencies w, we may use as 
the variational, or "objective" or "cost" function, for each distinct (positive) frequency component the 
spectral density associated with this expression: 

^S[a,a\3,r]{uj) = Jd^C [a* -{d^ + u^)a + a* -j + d-q] + ex., (242) 

whose variations, subject to the solenoidal gauge constraint and given spatial boundary conditions, lead 
to Euler-Lagrange equations consisting of the inhomogeneous Helmholtz equation for the vector potential 
and its complex conjugate at the frequency lo > 0. Choosing 

9 = 5^ [i (l - iSc) «x (l + iSc) + S ^C-^] , (243) 
or equivalently (at least asymptotically, as ^ ^ oo). 
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for some undetermined real parameters A = A(w) and Vq = Pq{uj), equatfon H242() becomes: 

£S[a,a*JJ*]{u;) = 2j£CRe[^*-i9^ + ''^)^ + ' + ^(^) 1"°"] («^; ^) " (^)] ■ (245) 

The chosen boundary term amounts to the addition of a Lagrange muhipher enforcing a specified out- 
going power spectral density in the radiation fields at the stationary points of -^S. 

Now, a natural approach to approximation, associated variously with the names of Raleigh, Ritz, and 
Galerkin in slightly different contexts, consists in simply replacing this infinite-dimensional variational 
problem with a finite-dimensional restriction or projection of it onto some more easily characterized 
space of possibilities. 

That is, we restrict the space of possible vector potentials to some parameterized family a((^; w; ex) for 
some finite-dimensional vector a of continuous parameters (over some specified domain), and thereby 
replace a variational problem involving functional derivatives of the action functional with respect to 
the vector potential with one involving ordinary derivatives of an action function with respect to the 
variational shape parameters determining the form of the radiation and the Lagrange multiplier enforcing 
a power normalization. That is, assuming the family of variational trial functions is implicitly constrained 
to satisfy the gauge constraint, the inhomogeneous Helmholtz PDE is replaced for each a; > with a 
system of algebraic equations, 

££5(«) = 0, (246a) 
^£S{c^) = 0, (246b) 

whose solution a — a[j]{u!;VQ), A — X[j]{uj;VQ) determines the variational approximation a{<^;uj;a.) to 
a(C;w). 

Now, following our earlier development, suppose we consider a family x(C; t^; ck) of trial vector fields for 
the general inhomogeneous problem which are explicitly constrained to be solenoidal and to satisfy the 
source-free Helmholtz equation, so the first term in the action density vanishes identically. Then the 
variational action spectral density becomes 

ijS{u;; a; A) - ^ Im {e^^{cj; a)\j{cj)) + X{u;) [|^y,M[x°"1(oo; u;) - 7'^(^)] . (247) 

Despite certain similarities, the Ritz-Galerkin-type variational principle associated with this action 
is distinct form the MPVP. The action principle involves finding constrained stationary points of 
Im. {e^ ±{<jj)\j {(^)) , while the MPVP involves finding the constrained maxima of T^^^-t.le-^ ±, j]{uj; a) = 
Re(£,t_L(w;a)|j(w)) . 

Although generically the "true" Euler-Lagrange solutions correspond to saddle-points, when the trial 
functions are restricted to a parameterized source-free family and when the power-normalization 
constraint is imposed via a Lagrange multiplier, the overall problem inherits convexity from the 
constraint alone, and the critical points turn out to be power-constrained local maxima or min- 
ima of lm{£~^±{uj;a)\j{Lu)) , which implies that the overall phase should be chosen such that 
Re {e-^±{uj; cy.)\j (uj)) — 0. In the MPVP case, the variational solution corresponds to constrained local 
maxima of — Re (£^_l((jj; q:)|j('^)) i so that Im (e-^ j^((jj; Q:)|j(aj)) — 0. Any attempt to employ directly 
the action-based variational principle with the source-free trial solutions thus encounters two related 
problems: the ostensible method converges to a solution with a global phase error of 7r/2; and as an 
immediate result, the calculated work performed on/by the sources vanishes, so the problem becomes 
degenerate, and we are left with no means to consistently choose the absolute power level. 

It is not difficult to trace the origin of this phase error. Naive application of Fourier transforms to the 
Green function yields the reciprocal-space, frequency-domain kernel 

G-\k,k-u)^5{k-k') ^ \ (248) 
oj'^ — ||«;|| 

but this neglects the requirement of causality demanded of the retarded Green function, wherein the 
response can appear only after the source is applied, implying that the Green function must be analytic 
for Imw > 0. We must be slightly more careful, and here treat the (temporal) Fourier transform as the 
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limiting case of a Laplace transform, where the initial conditions are pushed back into the arbitrarily 
remote past rather than future, before being forgotten altogether. The Green function (|248|l must be 
understood in the sense 



e^0+ ll'^ll — (^^ + 



(249) 



where the limit is taken only after G""' is convolved with the Fourier transform of the source. In the 
full Fourier representation, the vector potential and transverse source are then related by 



a{k; u) 



lim 



where 



and 



j^{k-u) = {\-kk')j{k-uj), 
j{k;uj) = J d^C j(C;c^)e-*'^ 



are the (scaled) Fourier transforms. For real oj and |1A;|| > \lu\ , we see that 

e a{k; us) 



while for llfcll < \lo\ 



arg 



arg 



e-j±{k;uj) 



e-a{k; uj) 
e-jj_{k;uj) 



0, 



(250) 



(251) 



(252a) 

(252b) 
(252c) 



(253) 



(254) 



In either case, the associated electric field for these spectral components is then in quadrature with 
the current, and no time-averaged work is performed or energy exchanged. However, on the singular 
resonant manifold where ||fc|| — \uj\ exactly, the magnitude of the Green function diverges in the limit as 
e ^ 0+, but the phase is such that 



arg 



e-a{k\u) = ||fc||) 
e-j_i{k;u; = \\k\\) 



(255) 



The relative phase between field and source then jumps discontinuously, to where the currents deliver 
time-averaged energy to the fields. Our purported action principle fails to capture this, while the 
MP VP, although operating entirely within this singular manifold of source- free solutions, relies on energy 
conservation, which then enforces the correct phase relation. This is not to imply that action-based 
variational principles are not also useful for radiation problems, only that they will not naturally work 
with source-free solutions, and generically will involve finding saddle-points, not extrema. 

In fact, most of the variational principles previously employed for FEL analysis, paraxial wave propaga- 
tion, or laser-plasma problems are based on precisely this type of generalized Ritz-Galerkin approxima- 
tions to the underlying Lagrangian dynamics, although some of the authors have erroneously claimed 
extremal instead of mere stationary character for their techniques. A Lagrangian-based variational prin- 
ciple is developed in for paraxial optical propagation in an inhoniogeneous gain medium, which 
is generalized in ,52] to include nonlinear self-focusing effects, but the authors erroneously claim ex- 
tremal, not just stationary, properties. This approach was further extended to include more general 
laser-plasma interaction terms in l53|. where the authors correctly state that their trial solutions de- 
rive from a stationary principle. In |49l | the authors state without proof an extremal principle for the 
fundamental FEL mode in the small-signal regime, which in fact appears to be perfectly correct, but 
higher modes will merely be stationary. These ideas are further generalized in |m . Issl Issl Is^ , where 



49 



the variational techniques are erroneously described as extremal principles. Similar trial- function-based 
variational principles are also used in |57l l58> i59j for which, correctly, only stationarity is claimed. 

A number of other authors have made extensive use of (correctly-stated) stationary action principles 
both for general wave- wave and wave-particle interactions in kinetic and fluid treatments of plasmas and 
particle beams "gT, |^ and more specifically for analysis of FELs 65]. However, rather than 

involving parameter-laden trial functions, these approaches employ the variational principles in order 
to provide economical descriptions and derivations of dynamical equations and conservation laws in the 
form of conventional ODEs or PDEs, and to effect in an efficient and transparent manner certain eikonal 
expansions or approximations, or oscillation-center/ponderomotive or other averaging procedures (66i.l67j| 
while preserving the Hamiltonian nature of the dynamics. 

Such confusion between stationary and extremal principles seems deeply embedded in physics literature 
and folklore. Linguistically, "optimizing," "maximizing," or "minimizing" simply sound better than 
"criticalizing" or "making stationary." Philosophically, optimality enjoys a certain teleological appeal 
which mere stationarity lacks. Historically, the buzzwords "least action" have been invoked so often 
in the discussion of dynamics that Hamilton's Principle and the Principle of Least Action have been 
mistakenly conflated, despite being quite distinct concepts. The Principle of Least Action was first 
articulated by Maupertuis and formalized by Euler, not Hamilton; is restricted to a smaller class of 
Lagrangians (no explicit time dependence); uses a different action (the abbreviated action integral) and 
a different set of constrained variations (the energy is fixed but temporal endpoints are not); and in fact 
it is itself misnamed, requiring only that the abbreviated action be stationary, not necessarily minimal, 
for the physical trajectory. Practically, both optimization principles and stationary principles share an 
insensitivity of the variational quantity to the trial functions, but the extremal case additionally provides 
an upper or lower bound. 

Numerically, more efficient and/or reliable computational techniques may be employed for optimization 
problems. Finding (local) maxima or minima of functions is much easier than finding roots of systems of 
equations. The "Alpine" analogy is often made: lost in a foggy terrain, the mountain-climber can reliably 
find a nearby peak (local maximum) by moving uphill, a nearby valley (local minimum) by wandering 
downhill, but has no guaranteed strategy for finding a mountain pass (saddle-point.) If the original 
action density is convex (at least in some sufficiently large region of function space), so that the true 
solution corresponds to an extremum of the action, then in the Raleigh-Ritz approximation we have a 
natural metric to measure both "closeness" to the true solution and "progress" in iteratively determining 
the approximate one - namely, a metric induced by the action itself. Without convexity, the true and 
approximate solutions may correspond to saddle-points, and we must introduce some arbitrary external 
metric to measure similarity or improvement. Even in the classic formulations which employ only linear 
expansions in basis functions, the extremal case leads to the solution of equations corresponding to a 
symmetric (or Hermitian), positive definite matrix, while in the merely stationary problem the matrix 
is in general Hermitian but not positive definite. 

In the case of electromagnetic radiation, the stationary points are generically saddles, because of the 
hyperbolic nature of the wave equation. Alternatively, one can see this by the well-known equivalence of 
electromagnetic mode dynamics to a collection of harmonic oscillators, for which the action functional is 
known to be minimal only for sufficiently short times intervals away from turning points along an orbit. 

Similar subtleties occur in the usual Raleigh-Ritz approximation for the stationary states in quantum 
mechanics. The energy-expectation value is always a local minimum for the ground state, and is always 
stationary for any excited state, but is only minimal for excited states if the variations are constrained 
to be orthogonal to all lower states. 

To summarize, the MPVP involves finding constrained extrema of quantities of the form: 



while Rumsey reaction-based variational principles would involve finding constrained stationary points 
of quantities of the form: 




(256) 




(257) 
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and Lagrangian action-based principles would involve finding constrained stationary points of quantities 
of the form: 



We have reviewed in some detail a Hilbert-space and operator-based approach to electromagnetic ra- 
diation, and have used this formalism to derive in some detail a maximum-power variational principle 
(MPVP) for spontaneous radiation from prescribed classical sources, first in the paraxial limit and then 
in a more general setting, adding to the large family of variational techniques for electromagnetic prob- 
lems in general, and undulator/wiggler radiation in particular. Mathematical details aside, at its most 
essential, the MPVP is really just a straightforward and rather obvious consequence of two simple and 
rather obvious constraints, together with another fundamental mathematical fact: the power in any one 
sourcc-frcc mode of the electromagnetic field may not exceed the total power in all the modes (i.e., Bcssel 
inequality); and the power radiated must be attributable to power delivered by the sources, even in the 
regime where we ignore back-action on the sources (i.e., conservation of energy); while information about 
the fields on some boimdary surface, needed to determine this radiated power, can be converted into in- 
formation about the derivatives of the fields in the interior (i.e., Gauss's law, one of the multidimensional 
generalizations of the Fundamental Theorem of Calculus). 

This approximation, or ones similar to it, has been frequently used, almost without comment or per- 
ceived need for further justification, in classical or quantum stimulated emission situations, where in the 
presence of gain we naturally expect to observe that mode which grows the fastest, but it is equally 
applicable in the spontaneous regime, because arguments along the lines of Einstein's derivation of the 
A and B coefficients lead to definite connections between spontaneous emission, stimulated emission, 
and stimulated absorption, even when the radiation is completely classical. However simple, even triv- 
ial, these observations are not without practical content or application to undulator and possibly other 
radiation problems. 
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